/**  ------ VALUE FUNCTION ITERATION FILE ------ */







// ADD NEEDED FILES //
#include "values/valueB.cpp"
#include "values/valueR.cpp"
#include "values/valueE.cpp"
#include "values/valueEexclu.cpp"
#include "values/valueEexcluhat.cpp"
#include "values/valueU_E.cpp"
#include "values/valueW_E.cpp"




/******************************/
/** VALUE FUNCTION ITERATION **/
/******************************/


// find value functions (Value function Iteration).
#if SEA == 1
void VFI(double *valueUiNE, double *saveUiNE, double *C_UiNE, double *searchUiNE_W, double *searchUiNE_J,
         double *valueUiE,  double *saveUiE,  double *C_UiE,  double *searchUiE_W,  double *searchUiE_J,
         double *valueUnNE, double *saveUnNE, double *C_UnNE, double *searchUnNE_W, double *searchUnNE_J,
         double *valueUnE,  double *saveUnE,  double *C_UnE,  double *searchUnE_W,  double *searchUnE_J,
         double *valueWWNE, double *saveWWNE, double *C_WWNE, double *searchWWNE_J,
         double *valueWWE,  double *saveWWE,  double *C_WWE,  double *searchWWE_J,
         double *valueJE,   double *saveJE,   double *C_JE,   double *searchJE_W,   double *collateralJE,
         double *valueJNE,  double *saveJNE,  double *C_JNE,  double *searchJNE_W,  double *collateralJNE, double *sloanJNE, double *rateJNE, double *defaultJNE,
         double *valueJEi,  double *saveJEi,  double *C_JEi,  double *searchJEi_W,  double *collateralJEi,
         double *valueJNEi, double *saveJNEi, double *C_JNEi, double *searchJNEi_W, double *collateralJNEi, double *sloanJNEi, double *rateJNEi, double *defaultJNEi)
{
#endif
    
#if NOPOLICY == 1
    void VFI(double *valueUiNE, double *saveUiNE,  double *C_UiNE, double *searchUiNE_W, double *searchUiNE_J,
             double *valueUiE,  double *saveUiE,   double *C_UiE,  double *searchUiE_W,  double *searchUiE_J,
             double *valueUnNE, double *saveUnNE,  double *C_UnNE, double *searchUnNE_W, double *searchUnNE_J,
             double *valueUnE,  double *saveUnE,   double *C_UnE,  double *searchUnE_W,  double *searchUnE_J,
             double *valueWWNE, double *saveWWNE,  double *C_WWNE, double *searchWWNE_J,
             double *valueWWE,  double *saveWWE,   double *C_WWE,  double *searchWWE_J,
             double *valueJE,   double *saveJE,    double *C_JE,   double *searchJE_W,   double *collateralJE,
             double *valueJNE,  double *saveJNE,   double *C_JNE,  double *searchJNE_W,  double *collateralJNE, double *sloanJNE, double *rateJNE, double *defaultJNE)
    {
#endif
        
        
        // GLOBAL //
        int tgrid, igrid, zgrid, zz, bb, znext, ygrid, igridCOH, ygridnext, tgridnext, zgridnext, iter, iexcluded, ilongrun, iinsured, icaseJ, icase, ixgridB, ixgridR;
        double C_expUW_JE = 0.0, C_expUW_JNE = 0.0, C_expUW_JEi = 0.0, C_expUW_JNEi = 0.0;
        double critereVF, tempE, COHB, COHR, valueBB, valueRR, valmax, bruteax, tmp_point, amaxAc, savemax, amaxsave, searchval_J, searchval_W, valfnmax, cohB, cohR, xgridB, dxgridB, xgridR, dxgridR, Bval, Rval, profitR, muENDO,RATEendo,collateral,C_tempE,C_BB,C_RR;
        
        // optimization values.
        double focE_grid, focE_gridm01, focE_min, focE_minp01, focE_max, focE_maxm01, RB_max;  // To find threshold at which he defaults //
        double valopt1, axopti1, rate_endo_opti; // 1st SOLUTION //
        double valopt2, axopti2, lowerbound, upperbound; // 2nd SOLUTION //
        double valopt3, axopti3, pointbb3_aft, pointbb3_bef; // 3rd SOLUTION //
        
        
        double *valueB, *valueR, *valueEhat;                  // values linked to default / repayment
        valueB          = (double *) calloc((ifulldimRB), sizeof(double));  // the dimension includes maxgridCOH
        valueR          = (double *) calloc((ifulldimRB), sizeof(double));  // the dimension includes maxgridCOH
        valueEhat       = (double *) calloc((ifulldimRB), sizeof(double));  // the dimension includes maxgridCOH
        double *C_B, *C_R, *C_Ehat;                  // values linked to default / repayment
        C_B          = (double *) calloc((ifulldimRB), sizeof(double));  // the dimension includes maxgridCOH
        C_R          = (double *) calloc((ifulldimRB), sizeof(double));  // the dimension includes maxgridCOH
        C_Ehat       = (double *) calloc((ifulldimRB), sizeof(double));  // the dimension includes maxgridCOH
        double *searchB, *searchR, *saveB, *saveR, *saveEhat, *searchEhat;
        searchB     = (double *) calloc((ifulldimRB), sizeof(double));
        searchR     = (double *) calloc((ifulldimRB), sizeof(double));
        saveB       = (double *) calloc((ifulldimRB), sizeof(double));
        saveR       = (double *) calloc((ifulldimRB), sizeof(double));
        saveEhat    = (double *) calloc((ifulldimRB), sizeof(double));
        searchEhat  = (double *) calloc((ifulldimRB), sizeof(double));
        double *valueJEz, *valueJNEz;
        valueJEz        = (double *) calloc((ifulldimEE), sizeof(double));
        valueJNEz       = (double *) calloc((ifulldimEE), sizeof(double));
        
        
        double *valueJNEtemp, *saveJNEtemp, *searchJNEtemp, *defaultJNEtemp, *valueJEtemp, *saveJEtemp, *searchJEtemp;
        valueJNEtemp    = (double *) calloc((maxfirmtype), sizeof(double));
        saveJNEtemp     = (double *) calloc((maxfirmtype), sizeof(double));
        searchJNEtemp   = (double *) calloc((maxfirmtype), sizeof(double));
        defaultJNEtemp  = (double *) calloc((maxfirmtype), sizeof(double));
        valueJEtemp     = (double *) calloc((maxfirmtype), sizeof(double));
        saveJEtemp      = (double *) calloc((maxfirmtype), sizeof(double));
        searchJEtemp    = (double *) calloc((maxfirmtype), sizeof(double));
        
#if SEA == 1
        double *valueBi, *valueRi, *valueEhati;
        valueBi          = (double *) calloc((ifulldimRB), sizeof(double));
        valueRi          = (double *) calloc((ifulldimRB), sizeof(double));
        valueEhati       = (double *) calloc((ifulldimRB), sizeof(double));
        double *C_Bi, *C_Ri, *C_Ehati;
        C_Bi          = (double *) calloc((ifulldimRB), sizeof(double));
        C_Ri          = (double *) calloc((ifulldimRB), sizeof(double));
        C_Ehati       = (double *) calloc((ifulldimRB), sizeof(double));
        double *searchBi, *searchRi, *saveBi, *saveRi, *saveEhati, *searchEhati;
        searchBi    = (double *) calloc((ifulldimRB), sizeof(double));
        searchRi    = (double *) calloc((ifulldimRB), sizeof(double));
        saveBi      = (double *) calloc((ifulldimRB), sizeof(double));
        saveRi      = (double *) calloc((ifulldimRB), sizeof(double));
        searchEhati = (double *) calloc((ifulldimRB), sizeof(double));
        saveEhati   = (double *) calloc((ifulldimRB), sizeof(double));
        double  *valueJEiz, *valueJNEiz;
        valueJEiz       =  (double *) calloc((ifulldimEE), sizeof(double));
        valueJNEiz      = (double *) calloc((ifulldimEE), sizeof(double));
#endif
        
        
        // expected value functions.
        double expUW_JE = 0.0, expUW_JNE = 0.0, expUW_JEi = 0.0, expUW_JNEi = 0.0;
        double *expE_UnE, *expE_UiE, *expE_WWE, *expE_WWNE_JNE, *expE_WWE_JE, *expE_UnNE_JNE, *expE_UnE_JE;
        double *expU_UnNE_JNE, *expU_UiNE_JNE, *expU_UnE_JE, *expU_UiE_JE, *expU_UnNE, *expU_UnE, *expU_UiNE, *expU_UiE, *expU_WWNE, *expU_WWE, *expU_WWNE_JNE, *expU_WWE_JE;
        double *expW_WWE_JE, *expW_WWNE_JNE, *expW_WWE, *expW_WWNE;
        expE_UnE        = (double *) calloc((ifulldimE), sizeof(double));
        expE_UiE        = (double *) calloc((ifulldimE), sizeof(double));
        expE_WWE        = (double *) calloc((ifulldimE), sizeof(double));
        expE_WWNE_JNE   = (double *) calloc((ifulldimE), sizeof(double));
        expE_WWE_JE     = (double *) calloc((ifulldimE), sizeof(double));
        expE_UnNE_JNE   = (double *) calloc((ifulldimE), sizeof(double));
        expE_UnE_JE     = (double *) calloc((ifulldimE), sizeof(double));
        
        expU_UnNE_JNE   = (double *) calloc((ifulldim), sizeof(double));
        expU_UiNE_JNE   = (double *) calloc((ifulldim), sizeof(double));
        expU_UnE_JE     = (double *) calloc((ifulldim), sizeof(double));
        expU_UiE_JE     = (double *) calloc((ifulldim), sizeof(double));
        expU_UnNE       = (double *) calloc((ifulldim), sizeof(double));
        expU_UnE        = (double *) calloc((ifulldim), sizeof(double));
        expU_UiNE       = (double *) calloc((ifulldim), sizeof(double));
        expU_UiE        = (double *) calloc((ifulldim), sizeof(double));
        expU_WWNE       = (double *) calloc((ifulldim), sizeof(double));
        expU_WWE        = (double *) calloc((ifulldim), sizeof(double));
        expU_WWNE_JNE   = (double *) calloc((ifulldim), sizeof(double));
        expU_WWE_JE     = (double *) calloc((ifulldim), sizeof(double));
        
        expW_WWE_JE     = (double *) calloc((ifulldimW), sizeof(double));
        expW_WWNE_JNE   = (double *) calloc((ifulldimW), sizeof(double));
        expW_WWE        = (double *) calloc((ifulldimW), sizeof(double));
        expW_WWNE       = (double *) calloc((ifulldimW), sizeof(double));
#if OPT_compute_C_flow == 1
        double *C_expE_UnE, *C_expE_UiE, *C_expE_WWE, *C_expE_WWNE_JNE, *C_expE_WWE_JE, *C_expE_UnNE_JNE, *C_expE_UnE_JE;
        double *C_expU_UnNE_JNE, *C_expU_UiNE_JNE, *C_expU_UnE_JE, *C_expU_UiE_JE, *C_expU_UnNE, *C_expU_UnE, *C_expU_UiNE, *C_expU_UiE, *C_expU_WWNE, *C_expU_WWE, *C_expU_WWNE_JNE, *C_expU_WWE_JE;
        double *C_expW_WWE_JE, *C_expW_WWNE_JNE, *C_expW_WWE, *C_expW_WWNE;
        C_expE_UnE        = (double *) calloc((ifulldimE), sizeof(double));
        C_expE_UiE        = (double *) calloc((ifulldimE), sizeof(double));
        C_expE_WWE        = (double *) calloc((ifulldimE), sizeof(double));
        C_expE_WWNE_JNE   = (double *) calloc((ifulldimE), sizeof(double));
        C_expE_WWE_JE     = (double *) calloc((ifulldimE), sizeof(double));
        C_expE_UnNE_JNE   = (double *) calloc((ifulldimE), sizeof(double));
        C_expE_UnE_JE     = (double *) calloc((ifulldimE), sizeof(double));
        
        C_expU_UnNE_JNE   = (double *) calloc((ifulldim), sizeof(double));
        C_expU_UiNE_JNE   = (double *) calloc((ifulldim), sizeof(double));
        C_expU_UnE_JE     = (double *) calloc((ifulldim), sizeof(double));
        C_expU_UiE_JE     = (double *) calloc((ifulldim), sizeof(double));
        C_expU_UnNE       = (double *) calloc((ifulldim), sizeof(double));
        C_expU_UnE        = (double *) calloc((ifulldim), sizeof(double));
        C_expU_UiNE       = (double *) calloc((ifulldim), sizeof(double));
        C_expU_UiE        = (double *) calloc((ifulldim), sizeof(double));
        C_expU_WWNE       = (double *) calloc((ifulldim), sizeof(double));
        C_expU_WWE        = (double *) calloc((ifulldim), sizeof(double));
        C_expU_WWNE_JNE   = (double *) calloc((ifulldim), sizeof(double));
        C_expU_WWE_JE     = (double *) calloc((ifulldim), sizeof(double));
        
        C_expW_WWE_JE     = (double *) calloc((ifulldimW), sizeof(double));
        C_expW_WWNE_JNE   = (double *) calloc((ifulldimW), sizeof(double));
        C_expW_WWE        = (double *) calloc((ifulldimW), sizeof(double));
        C_expW_WWNE       = (double *) calloc((ifulldimW), sizeof(double));
#endif
        
#if SEA == 1
        double *expE_WWNE_JNEi, *expE_WWE_JEi, *expE_UiE_JEi, *expE_UiNE_JNEi;
        double *expU_UiNE_JNEi, *expU_UiE_JEi, *expU_WWNE_JNEi, *expU_WWE_JEi;
        expE_WWNE_JNEi      = (double *) calloc((ifulldimE), sizeof(double));
        expE_WWE_JEi        = (double *) calloc((ifulldimE), sizeof(double));
        expE_UiE_JEi        = (double *) calloc((ifulldimE), sizeof(double));
        expE_UiNE_JNEi      = (double *) calloc((ifulldimE), sizeof(double));
        expU_UiNE_JNEi      = (double *) calloc((ifulldim), sizeof(double));
        expU_UiE_JEi        = (double *) calloc((ifulldim), sizeof(double));
        expU_WWNE_JNEi      = (double *) calloc((ifulldim), sizeof(double));
        expU_WWE_JEi        = (double *) calloc((ifulldim), sizeof(double));
#if OPT_compute_C_flow == 1
        double *C_expE_WWNE_JNEi, *C_expE_WWE_JEi, *C_expE_UiE_JEi, *C_expE_UiNE_JNEi;
        double *C_expU_UiNE_JNEi, *C_expU_UiE_JEi, *C_expU_WWNE_JNEi, *C_expU_WWE_JEi;
        C_expE_WWNE_JNEi      = (double *) calloc((ifulldimE), sizeof(double));
        C_expE_WWE_JEi        = (double *) calloc((ifulldimE), sizeof(double));
        C_expE_UiE_JEi        = (double *) calloc((ifulldimE), sizeof(double));
        C_expE_UiNE_JNEi      = (double *) calloc((ifulldimE), sizeof(double));
        C_expU_UiNE_JNEi      = (double *) calloc((ifulldim), sizeof(double));
        C_expU_UiE_JEi        = (double *) calloc((ifulldim), sizeof(double));
        C_expU_WWNE_JNEi      = (double *) calloc((ifulldim), sizeof(double));
        C_expU_WWE_JEi        = (double *) calloc((ifulldim), sizeof(double));
#endif
#endif
        
        
        // STRUCTURE //
        struct paramsgolden params;         // For solving worker (C-S-S problem), entrepreneur (k problem) and unemployed (C-S problem)
        struct paramsgoldenRB paramsRB;     // For solving bankruptcy (C-S-S problem) and Repayment (C-S-S problem)
        
        double howard_off = 0.0;
        critereVF=1.0;  iter=0;
#if indexPM == 0    // if transition
        while ((critereVF > epsilonValue)  && (iter < maxiterVF)){
#endif
            
            critereVF = 0.0;
            
            /***********************************************************************/
            /** FIRST, COMPUTE THE EXPECTATIONS **/
            /***********************************************************************/
#if OMP == 1
#pragma omp parallel
            {
#pragma omp for private(igrid,tgrid,zgrid,ygrid,tgridnext,ygridnext,znext,expUW_JNE,expUW_JE,expUW_JEi,expUW_JNEi, C_expUW_JNE,C_expUW_JE,C_expUW_JEi,C_expUW_JNEi)
#endif
                for(igrid=0;igrid<maxgrid;igrid++){
                    for(tgrid=0;tgrid<maxtheta;tgrid++){
                        for(zgrid=0;zgrid<maxfirmtype;zgrid++) {
                            
                            // ********  Expectations of the Entrepreneurs.
                            expE_UnE[inxE(igrid,tgrid,zgrid)]       = 0.0;  expE_UiE[inxE(igrid,tgrid,zgrid)]     = 0.0;    expE_WWE[inxE(igrid,tgrid,zgrid)]       = 0.0;
                            expE_WWNE_JNE[inxE(igrid,tgrid,zgrid)]  = 0.0;  expE_WWE_JE[inxE(igrid,tgrid,zgrid)]  = 0.0;    expE_UnNE_JNE[inxE(igrid,tgrid,zgrid)]  = 0.0;
                            expE_UnE_JE[inxE(igrid,tgrid,zgrid)]    = 0.0;
                            
#if OPT_compute_C_flow == 1
                            C_expE_UnE[inxE(igrid,tgrid,zgrid)]       = 0.0;  C_expE_UiE[inxE(igrid,tgrid,zgrid)]     = 0.0;    C_expE_WWE[inxE(igrid,tgrid,zgrid)]       = 0.0;
                            C_expE_WWNE_JNE[inxE(igrid,tgrid,zgrid)]  = 0.0;  C_expE_WWE_JE[inxE(igrid,tgrid,zgrid)]  = 0.0;    C_expE_UnNE_JNE[inxE(igrid,tgrid,zgrid)]  = 0.0;
                            C_expE_UnE_JE[inxE(igrid,tgrid,zgrid)]    = 0.0;
#endif
#if SEA == 1
                            expE_WWNE_JNEi[inxE(igrid,tgrid,zgrid)] = 0.0;  expE_WWE_JEi[inxE(igrid,tgrid,zgrid)] = 0.0;    expE_UiE_JEi[inxE(igrid,tgrid,zgrid)]   = 0.0; expE_UiNE_JNEi[inxE(igrid,tgrid,zgrid)] = 0.0;
#if OPT_compute_C_flow == 1
                            C_expE_WWNE_JNEi[inxE(igrid,tgrid,zgrid)] = 0.0;  C_expE_WWE_JEi[inxE(igrid,tgrid,zgrid)] = 0.0;    C_expE_UiE_JEi[inxE(igrid,tgrid,zgrid)]   = 0.0; C_expE_UiNE_JNEi[inxE(igrid,tgrid,zgrid)] = 0.0;
#endif
#endif
                            
                            for(tgridnext=0;tgridnext<maxtheta;tgridnext++){
                                for(ygridnext = 0.0; ygridnext<maxyprod; ygridnext++) {
                                    
                                    expE_WWE[inxE(igrid,tgrid,zgrid)]   += invmyprod[ygridnext]*mprod[tgrid][tgridnext]*valueWWE[inxW(igrid,tgridnext,ygridnext)];
                                    
                                    expE_WWNE_JNE[inxE(igrid,tgrid,zgrid)] += invmyprod[ygridnext]*mprod[tgrid][tgridnext]*max(valueWWNE[inxW(igrid,tgridnext,ygridnext)],valueJNE[inxE(igrid,tgridnext,zgrid)]);
                                    expE_WWE_JE[inxE(igrid,tgrid,zgrid)]   += invmyprod[ygridnext]*mprod[tgrid][tgridnext]*max(valueWWE[inxW(igrid,tgridnext,ygridnext)],valueJE[inxE(igrid,tgridnext,zgrid)]);
#if OPT_compute_C_flow == 1
                                    C_expE_WWNE_JNE[inxE(igrid,tgrid,zgrid)] += invmyprod[ygridnext]*mprod[tgrid][tgridnext]*((valueWWNE[inxW(igrid,tgridnext,ygridnext)] >= valueJNE[inxE(igrid,tgridnext,zgrid)])?(C_WWNE[inxW(igrid,tgridnext,ygridnext)]):(C_JNE[inxE(igrid,tgridnext,zgrid)]));
                                    C_expE_WWE_JE[inxE(igrid,tgrid,zgrid)]   += invmyprod[ygridnext]*mprod[tgrid][tgridnext]*((valueWWE[inxW(igrid,tgridnext,ygridnext)] >= valueJE[inxE(igrid,tgridnext,zgrid)])?(C_WWE[inxW(igrid,tgridnext,ygridnext)]):(C_JE[inxE(igrid,tgridnext,zgrid)]));
                                    C_expE_WWE[inxE(igrid,tgrid,zgrid)] += invmyprod[ygridnext]*mprod[tgrid][tgridnext]*C_WWE[inxW(igrid,tgridnext,ygridnext)];
#endif
                                    
#if SEA == 1
                                    expE_WWNE_JNEi[inxE(igrid,tgrid,zgrid)] += invmyprod[ygridnext]*mprod[tgrid][tgridnext]*max(valueWWNE[inxW(igrid,tgridnext,ygridnext)],valueJNEi[inxE(igrid,tgridnext,zgrid)]);
                                    expE_WWE_JEi[inxE(igrid,tgrid,zgrid)]   += invmyprod[ygridnext]*mprod[tgrid][tgridnext]*max(valueWWE[inxW(igrid,tgridnext,ygridnext)],valueJEi[inxE(igrid,tgridnext,zgrid)]);
#if OPT_compute_C_flow == 1
                                    C_expE_WWNE_JNEi[inxE(igrid,tgrid,zgrid)] += invmyprod[ygridnext]*mprod[tgrid][tgridnext]*((valueWWNE[inxW(igrid,tgridnext,ygridnext)] >= valueJNEi[inxE(igrid,tgridnext,zgrid)])?(C_WWNE[inxW(igrid,tgridnext,ygridnext)]):(C_JNEi[inxE(igrid,tgridnext,zgrid)]));
                                    C_expE_WWE_JEi[inxE(igrid,tgrid,zgrid)]   += invmyprod[ygridnext]*mprod[tgrid][tgridnext]*((valueWWE[inxW(igrid,tgridnext,ygridnext)] >= valueJEi[inxE(igrid,tgridnext,zgrid)])?(C_WWE[inxW(igrid,tgridnext,ygridnext)]):(C_JEi[inxE(igrid,tgridnext,zgrid)]));
#endif
#endif
                                    
                                }
                                
                                expE_UnE[inxE(igrid,tgrid,zgrid)]       += mprod[tgrid][tgridnext]*valueUnE[inx(igrid,tgridnext)];
                                expE_UiE[inxE(igrid,tgrid,zgrid)]       += mprod[tgrid][tgridnext]*valueUiE[inx(igrid,tgridnext)];
                                
                                
                                expE_UnNE_JNE[inxE(igrid,tgrid,zgrid)]  += mprod[tgrid][tgridnext]*max(valueUnNE[inx(igrid,tgridnext)],valueJNE[inxE(igrid,tgridnext,zgrid)]);
                                expE_UnE_JE[inxE(igrid,tgrid,zgrid)]    += mprod[tgrid][tgridnext]*max(valueUnE[inx(igrid,tgridnext)],valueJE[inxE(igrid,tgridnext,zgrid)]);
                                
#if OPT_compute_C_flow == 1
                                C_expE_UnNE_JNE[inxE(igrid,tgrid,zgrid)]  += mprod[tgrid][tgridnext]*((valueUnNE[inx(igrid,tgridnext)] >= valueJNE[inxE(igrid,tgridnext,zgrid)])?(C_UnNE[inx(igrid,tgridnext)]):(C_JNE[inxE(igrid,tgridnext,zgrid)]));
                                C_expE_UnE_JE[inxE(igrid,tgrid,zgrid)]    += mprod[tgrid][tgridnext]*((valueUnE[inx(igrid,tgridnext)] >= valueJE[inxE(igrid,tgridnext,zgrid)])?(C_UnE[inx(igrid,tgridnext)]):(C_JE[inxE(igrid,tgridnext,zgrid)]));
                                C_expE_UnE[inxE(igrid,tgrid,zgrid)]     += mprod[tgrid][tgridnext]*C_UnE[inx(igrid,tgridnext)];
                                C_expE_UiE[inxE(igrid,tgrid,zgrid)]     += mprod[tgrid][tgridnext]*C_UiE[inx(igrid,tgridnext)];
#endif
                                
                                // expectation over t, knowing igrid and zgrid (realized during the period).
#if SEA == 1
                                expE_UiE_JEi[inxE(igrid,tgrid,zgrid)]   += mprod[tgrid][tgridnext]*max(valueUiE[inx(igrid,tgridnext)],valueJEi[inxE(igrid,tgridnext,zgrid)]);
                                expE_UiNE_JNEi[inxE(igrid,tgrid,zgrid)] += mprod[tgrid][tgridnext]*max(valueUiNE[inx(igrid,tgridnext)],valueJNEi[inxE(igrid,tgridnext,zgrid)]);
#if OPT_compute_C_flow == 1
                                C_expE_UiE_JEi[inxE(igrid,tgrid,zgrid)]   += mprod[tgrid][tgridnext]*((valueUiE[inx(igrid,tgridnext)] >= valueJEi[inxE(igrid,tgridnext,zgrid)])?(C_UiE[inx(igrid,tgridnext)]):(C_JEi[inxE(igrid,tgridnext,zgrid)]));
                                C_expE_UiNE_JNEi[inxE(igrid,tgrid,zgrid)] += mprod[tgrid][tgridnext]*((valueUiNE[inx(igrid,tgridnext)] >= valueJNEi[inxE(igrid,tgridnext,zgrid)])?(C_UiNE[inx(igrid,tgridnext)]):(C_JNEi[inxE(igrid,tgridnext,zgrid)]));
#endif
#endif
                            }
                        } // end zgrid
                        
                        
                        // ********  Expectations of the Unemployed
                        expU_UnNE_JNE[inx(igrid,tgrid)] = 0.0; expU_UiNE_JNE[inx(igrid,tgrid)] = 0.0; expU_UnE_JE[inx(igrid,tgrid)] = 0.0; expU_UiE_JE[inx(igrid,tgrid)]   = 0.0;
                        expU_UnNE[inx(igrid,tgrid)] = 0.0; expU_UnE[inx(igrid,tgrid)]  = 0.0; expU_UiNE[inx(igrid,tgrid)] = 0.0; expU_UiE[inx(igrid,tgrid)]  = 0.0;
                        expU_WWNE[inx(igrid,tgrid)] = 0.0; expU_WWE[inx(igrid,tgrid)]  = 0.0;
                        expU_WWNE_JNE[inx(igrid,tgrid)] = 0.0; expU_WWE_JE[inx(igrid,tgrid)]  = 0.0;
                        
#if OPT_compute_C_flow == 1
                        C_expU_UnNE_JNE[inx(igrid,tgrid)] = 0.0; C_expU_UiNE_JNE[inx(igrid,tgrid)] = 0.0; C_expU_UnE_JE[inx(igrid,tgrid)] = 0.0; C_expU_UiE_JE[inx(igrid,tgrid)]   = 0.0;
                        C_expU_UnNE[inx(igrid,tgrid)] = 0.0; C_expU_UnE[inx(igrid,tgrid)]  = 0.0; C_expU_UiNE[inx(igrid,tgrid)] = 0.0; C_expU_UiE[inx(igrid,tgrid)]  = 0.0;
                        C_expU_WWNE[inx(igrid,tgrid)] = 0.0; C_expU_WWE[inx(igrid,tgrid)]  = 0.0;
                        C_expU_WWNE_JNE[inx(igrid,tgrid)] = 0.0; C_expU_WWE_JE[inx(igrid,tgrid)]  = 0.0;
#endif
                        
#if SEA == 1
                        expU_UiNE_JNEi[inx(igrid,tgrid)] = 0.0; expU_UiE_JEi[inx(igrid,tgrid)]   = 0.0;
                        expU_WWNE_JNEi[inx(igrid,tgrid)] = 0.0; expU_WWE_JEi[inx(igrid,tgrid)]   = 0.0;
#if OPT_compute_C_flow == 1
                        C_expU_UiNE_JNEi[inx(igrid,tgrid)] = 0.0; C_expU_UiE_JEi[inx(igrid,tgrid)]   = 0.0;
                        C_expU_WWNE_JNEi[inx(igrid,tgrid)] = 0.0; C_expU_WWE_JEi[inx(igrid,tgrid)]   = 0.0;
#endif
#endif
                        
                        for(tgridnext=0;tgridnext<maxtheta;tgridnext++){
                            
                            expUW_JE = 0.0; expUW_JNE = 0.0;
                            C_expUW_JE = 0.0; C_expUW_JNE = 0.0;
#if SEA == 1
                            expUW_JEi = 0.0; expUW_JNEi = 0.0;
                            C_expUW_JEi = 0.0; C_expUW_JNEi = 0.0;
#endif
                            
                            for(znext=0;znext<maxfirmtype;znext++) {
                                expUW_JE  += mzinv[znext]*valueJE[inxE(igrid,tgridnext,znext)];
                                expUW_JNE += mzinv[znext]*valueJNE[inxE(igrid,tgridnext,znext)];
#if OPT_compute_C_flow == 1
                                C_expUW_JE  += mzinv[znext]*C_JE[inxE(igrid,tgridnext,znext)];
                                C_expUW_JNE += mzinv[znext]*C_JNE[inxE(igrid,tgridnext,znext)];
#endif
                                
#if SEA == 1
                                expUW_JEi  += mzinv[znext]*valueJEi[inxE(igrid,tgridnext,znext)];
                                expUW_JNEi += mzinv[znext]*valueJNEi[inxE(igrid,tgridnext,znext)];
#if OPT_compute_C_flow == 1
                                C_expUW_JEi  += mzinv[znext]*C_JEi[inxE(igrid,tgridnext,znext)];
                                C_expUW_JNEi += mzinv[znext]*C_JNEi[inxE(igrid,tgridnext,znext)];
#endif
#endif
                            }
                            
                            // ********  Expectations of the Unemployed
                            expU_UnNE_JNE[inx(igrid,tgrid)] += mprod[tgrid][tgridnext]*max(valueUnNE[inx(igrid,tgridnext)],expUW_JNE);
                            expU_UiNE_JNE[inx(igrid,tgrid)] += mprod[tgrid][tgridnext]*max(valueUiNE[inx(igrid,tgridnext)],expUW_JNE);
                            expU_UnE_JE[inx(igrid,tgrid)]   += mprod[tgrid][tgridnext]*max(valueUnE[inx(igrid,tgridnext)],expUW_JE);
                            expU_UiE_JE[inx(igrid,tgrid)]   += mprod[tgrid][tgridnext]*max(valueUiE[inx(igrid,tgridnext)],expUW_JE);
                            
                            
                            expU_UnNE[inx(igrid,tgrid)] += mprod[tgrid][tgridnext]*valueUnNE[inx(igrid,tgridnext)];
                            expU_UnE[inx(igrid,tgrid)]  += mprod[tgrid][tgridnext]*valueUnE[inx(igrid,tgridnext)];
                            expU_UiNE[inx(igrid,tgrid)] += mprod[tgrid][tgridnext]*valueUiNE[inx(igrid,tgridnext)];
                            expU_UiE[inx(igrid,tgrid)]  += mprod[tgrid][tgridnext]*valueUiE[inx(igrid,tgridnext)];
                            
#if OPT_compute_C_flow == 1
                            C_expU_UnNE_JNE[inx(igrid,tgrid)] += mprod[tgrid][tgridnext]*((valueUnNE[inx(igrid,tgridnext)] >= expUW_JNE)?(C_UnNE[inx(igrid,tgridnext)]):(C_expUW_JNE));
                            C_expU_UiNE_JNE[inx(igrid,tgrid)] += mprod[tgrid][tgridnext]*((valueUiNE[inx(igrid,tgridnext)] >= expUW_JNE)?(C_UiNE[inx(igrid,tgridnext)]):(C_expUW_JNE));
                            C_expU_UnE_JE[inx(igrid,tgrid)]   += mprod[tgrid][tgridnext]*((valueUnE[inx(igrid,tgridnext)] >= expUW_JE)?(C_UnE[inx(igrid,tgridnext)]):(C_expUW_JE));
                            C_expU_UiE_JE[inx(igrid,tgrid)]   += mprod[tgrid][tgridnext]*((valueUiE[inx(igrid,tgridnext)] >= expUW_JE)?(C_UiE[inx(igrid,tgridnext)]):(C_expUW_JE));
                            
                            C_expU_UnNE[inx(igrid,tgrid)] += mprod[tgrid][tgridnext]*C_UnNE[inx(igrid,tgridnext)];
                            C_expU_UnE[inx(igrid,tgrid)]  += mprod[tgrid][tgridnext]*C_UnE[inx(igrid,tgridnext)];
                            C_expU_UiNE[inx(igrid,tgrid)] += mprod[tgrid][tgridnext]*C_UiNE[inx(igrid,tgridnext)];
                            C_expU_UiE[inx(igrid,tgrid)]  += mprod[tgrid][tgridnext]*C_UiE[inx(igrid,tgridnext)];
#endif
                            
#if SEA == 1
                            expU_UiNE_JNEi[inx(igrid,tgrid)] += mprod[tgrid][tgridnext]*max(valueUiNE[inx(igrid,tgridnext)],expUW_JNEi);
                            expU_UiE_JEi[inx(igrid,tgrid)]   += mprod[tgrid][tgridnext]*max(valueUiE[inx(igrid,tgridnext)],expUW_JEi);
#if OPT_compute_C_flow == 1
                            C_expU_UiNE_JNEi[inx(igrid,tgrid)] += mprod[tgrid][tgridnext]*((valueUiNE[inx(igrid,tgridnext)] >= expUW_JNEi)?(C_UiNE[inx(igrid,tgridnext)]):(C_expUW_JNEi));
                            C_expU_UiE_JEi[inx(igrid,tgrid)]   += mprod[tgrid][tgridnext]*((valueUiE[inx(igrid,tgridnext)] >= expUW_JEi)?(C_UiE[inx(igrid,tgridnext)]):(C_expUW_JEi));
#endif
#endif
                            
                            // note the shocks is known at the moment of the choice.
                            for(ygridnext = 0.0; ygridnext<maxyprod; ygridnext++) {
                                expU_WWNE[inx(igrid,tgrid)] += invmyprod[ygridnext]*mprod[tgrid][tgridnext]*valueWWNE[inxW(igrid,tgridnext,ygridnext)];
                                expU_WWE[inx(igrid,tgrid)]  += invmyprod[ygridnext]*mprod[tgrid][tgridnext]*valueWWE[inxW(igrid,tgridnext,ygridnext)];
                                
                                expU_WWNE_JNE[inx(igrid,tgrid)] += invmyprod[ygridnext]*mprod[tgrid][tgridnext]*max(valueWWNE[inxW(igrid,tgridnext,ygridnext)],expUW_JNE);
                                expU_WWE_JE[inx(igrid,tgrid)]   += invmyprod[ygridnext]*mprod[tgrid][tgridnext]*max(valueWWE[inxW(igrid,tgridnext,ygridnext)],expUW_JE);
                                
#if OPT_compute_C_flow == 1
                                C_expU_WWNE[inx(igrid,tgrid)] += invmyprod[ygridnext]*mprod[tgrid][tgridnext]*C_WWNE[inxW(igrid,tgridnext,ygridnext)];
                                C_expU_WWE[inx(igrid,tgrid)]  += invmyprod[ygridnext]*mprod[tgrid][tgridnext]*C_WWE[inxW(igrid,tgridnext,ygridnext)];
                                
                                C_expU_WWNE_JNE[inx(igrid,tgrid)] += invmyprod[ygridnext]*mprod[tgrid][tgridnext]*((valueWWNE[inxW(igrid,tgridnext,ygridnext)] >= expUW_JNE)?(C_WWNE[inxW(igrid,tgridnext,ygridnext)]):(C_expUW_JNE));
                                C_expU_WWE_JE[inx(igrid,tgrid)]   += invmyprod[ygridnext]*mprod[tgrid][tgridnext]*((valueWWE[inxW(igrid,tgridnext,ygridnext)] >= expUW_JE)?(C_WWE[inxW(igrid,tgridnext,ygridnext)]):(C_expUW_JE));
#endif
                                
                                
#if SEA == 1
                                expU_WWNE_JNEi[inx(igrid,tgrid)] += invmyprod[ygridnext]*mprod[tgrid][tgridnext]*max(valueWWNE[inxW(igrid,tgridnext,ygridnext)],expUW_JNEi);
                                expU_WWE_JEi[inx(igrid,tgrid)]   += invmyprod[ygridnext]*mprod[tgrid][tgridnext]*max(valueWWE[inxW(igrid,tgridnext,ygridnext)],expUW_JEi);
#if OPT_compute_C_flow == 1
                                C_expU_WWNE_JNEi[inx(igrid,tgrid)] += invmyprod[ygridnext]*mprod[tgrid][tgridnext]*((valueWWNE[inxW(igrid,tgridnext,ygridnext)] >= expUW_JNEi)?(C_WWNE[inxW(igrid,tgridnext,ygridnext)]):(C_expUW_JNEi));
                                C_expU_WWE_JEi[inx(igrid,tgrid)]   += invmyprod[ygridnext]*mprod[tgrid][tgridnext]*((valueWWE[inxW(igrid,tgridnext,ygridnext)] >= expUW_JEi)?(C_WWE[inxW(igrid,tgridnext,ygridnext)]):(C_expUW_JEi));
#endif
#endif
                            }
                        }
                        
                        
                        // ********  Expectations of the Worker
                        for(ygrid = 0.0; ygrid<maxyprod; ygrid++) {
                            
                            expW_WWE_JE[inxW(igrid,tgrid,ygrid)]     = 0.0; expW_WWNE_JNE[inxW(igrid,tgrid,ygrid)]    = 0.0;
                            expW_WWE[inxW(igrid,tgrid,ygrid)]        = 0.0; expW_WWNE[inxW(igrid,tgrid,ygrid)]       = 0.0;
#if OPT_compute_C_flow == 1
                            C_expW_WWE_JE[inxW(igrid,tgrid,ygrid)]     = 0.0; C_expW_WWNE_JNE[inxW(igrid,tgrid,ygrid)]    = 0.0;
                            C_expW_WWE[inxW(igrid,tgrid,ygrid)]        = 0.0; C_expW_WWNE[inxW(igrid,tgrid,ygrid)]       = 0.0;
#endif
                            for(tgridnext=0;tgridnext<maxtheta;tgridnext++){
                                
                                expUW_JE = 0.0; expUW_JNE = 0.0; C_expUW_JE = 0.0; C_expUW_JNE = 0.0;
                                
                                for(znext=0;znext<maxfirmtype;znext++) {
                                    expUW_JE  += mzinv[znext]*valueJE[inxE(igrid,tgridnext,znext)];
                                    expUW_JNE += mzinv[znext]*valueJNE[inxE(igrid,tgridnext,znext)];
#if OPT_compute_C_flow == 1
                                    C_expUW_JE  += mzinv[znext]*C_JE[inxE(igrid,tgridnext,znext)];
                                    C_expUW_JNE += mzinv[znext]*C_JNE[inxE(igrid,tgridnext,znext)];
#endif
                                }
                                
                                for(ygridnext = 0.0; ygridnext<maxyprod; ygridnext++) {
                                    expW_WWE_JE[inxW(igrid,tgrid,ygrid)]     += myprod[ygrid][ygridnext]*mprod[tgrid][tgridnext]*max(valueWWE[inxW(igrid,tgridnext,ygridnext)],expUW_JE);
                                    expW_WWNE_JNE[inxW(igrid,tgrid,ygrid)]   += myprod[ygrid][ygridnext]*mprod[tgrid][tgridnext]*max(valueWWNE[inxW(igrid,tgridnext,ygridnext)],expUW_JNE);
                                    
                                    expW_WWE[inxW(igrid,tgrid,ygrid)]        += myprod[ygrid][ygridnext]*mprod[tgrid][tgridnext]*valueWWE[inxW(igrid,tgridnext,ygridnext)];
                                    expW_WWNE[inxW(igrid,tgrid,ygrid)]       += myprod[ygrid][ygridnext]*mprod[tgrid][tgridnext]*valueWWNE[inxW(igrid,tgridnext,ygridnext)];
                                    
#if OPT_compute_C_flow == 1
                                    C_expW_WWE_JE[inxW(igrid,tgrid,ygrid)]     += myprod[ygrid][ygridnext]*mprod[tgrid][tgridnext]*((valueWWE[inxW(igrid,tgridnext,ygridnext)] >= expUW_JE)?(C_WWE[inxW(igrid,tgridnext,ygridnext)]):(C_expUW_JE));
                                    C_expW_WWNE_JNE[inxW(igrid,tgrid,ygrid)]   += myprod[ygrid][ygridnext]*mprod[tgrid][tgridnext]*((valueWWNE[inxW(igrid,tgridnext,ygridnext)] >= expUW_JNE)?(C_WWNE[inxW(igrid,tgridnext,ygridnext)]):(C_expUW_JNE));
                                    C_expW_WWE[inxW(igrid,tgrid,ygrid)]        += myprod[ygrid][ygridnext]*mprod[tgrid][tgridnext]*C_WWE[inxW(igrid,tgridnext,ygridnext)];
                                    C_expW_WWNE[inxW(igrid,tgrid,ygrid)]       += myprod[ygrid][ygridnext]*mprod[tgrid][tgridnext]*C_WWNE[inxW(igrid,tgridnext,ygridnext)];
#endif
                                    
                                } // end ygridnext
                            } // end tgridnext
                        } // end ygrid
                    } // end tgrid
                } // end igrid
#if OMP == 1
            } // end pragma
#endif
            
            
            /***********************************************************************/
            /** SECOND, solve the consumption-saving problem of entrepreneur
             /***********************************************************************/
#if OMP == 1
#pragma omp parallel
            {
#pragma omp for private(paramsRB, tgrid, igrid, zgrid, zz, bb, znext, ygrid, igridCOH, ygridnext, tgridnext, zgridnext, iter, iexcluded, ilongrun, iinsured, icaseJ, icase, ixgridB, ixgridR, valmax, bruteax, amaxAc, savemax, amaxsave, searchval_J, searchval_W, valfnmax, cohB, cohR, xgridB, dxgridB, xgridR, dxgridR, Bval, Rval, focE_grid, focE_gridm01, focE_min, focE_minp01, focE_max, focE_maxm01, RB_max,valopt1, axopti1, rate_endo_opti,valopt2, axopti2, lowerbound, upperbound,valopt3, axopti3, pointbb3_aft, pointbb3_bef)
#endif
                // -- BLOCK OK ----
                for(igridCOH=0; igridCOH<maxgridCOH; igridCOH++){       // solved on a grid of COH.
                    for(tgrid=0;tgrid<maxtheta;tgrid++){
                        for(zz=0;zz<maxfirmtype;zz++){
                            
                            
                            // PARAMS for structures //
                            paramsRB.JinsuredST = 0;
                            paramsRB.igridCOHST = igridCOH; paramsRB.tgridST = tgrid; paramsRB.zgridST = zz;
                            
                            // pass on the expected values.
                            paramsRB.expE_WWNE_JNE_ST   = expE_WWNE_JNE;
                            paramsRB.expE_WWE_JE_ST     = expE_WWE_JE;
                            paramsRB.expE_UnNE_JNE_ST   = expE_UnNE_JNE;
                            paramsRB.expE_UnE_JE_ST     = expE_UnE_JE;
                            paramsRB.expE_UnE_ST        = expE_UnE;
                            paramsRB.expE_WWE_ST        = expE_WWE;
#if OPT_compute_C_flow == 1
                            paramsRB.C_expE_WWNE_JNE_ST   = C_expE_WWNE_JNE;
                            paramsRB.C_expE_WWE_JE_ST     = C_expE_WWE_JE;
                            paramsRB.C_expE_UnNE_JNE_ST   = C_expE_UnNE_JNE;
                            paramsRB.C_expE_UnE_JE_ST     = C_expE_UnE_JE;
                            paramsRB.C_expE_UnE_ST        = C_expE_UnE;
                            paramsRB.C_expE_WWE_ST        = C_expE_WWE;
#endif
                            
                            /** ----  Solving Excluded HAT ---- **/
                            icase = 100; 
                            savemax = min(Gridmax-0.000000001,findconsRB((0.0),gridCOH[igridCOH])-0.000000001);
                            if(igridCOH == 0) {
                                valfnmax = focEhat(0.0,&paramsRB); 
                                amaxsave = Gridmin; searchval_W = paramsRB.searchWoutST;
                                icase = 1.0;
                            } else {
                                if(focEhat(Gridmin,&paramsRB) < focEhat(Gridmin+0.0000001, &paramsRB)) {
                                    valfnmax = focEhat(0.0,&paramsRB);  
                                    searchval_W = paramsRB.searchWoutST;  amaxsave = 0.0;
                                    icase = 1.0;
                                }else{//cas ou le max est a droite de la borne inf
                                    
                                    if(focEhat(savemax,&paramsRB) < focEhat((savemax-0.0000001), &paramsRB)) {
                                        valfnmax = focEhat(savemax,&paramsRB); 
                                        searchval_W = paramsRB.searchWoutST; amaxsave = savemax;
                                        icase = 3.0;
                                    }
                                    
                                    if(icase >= 100) {  // interior case.
                                        valfnmax = mygolden(GridminCOH,(GridminCOH+0.0000001),savemax,TOL,amaxsave,&paramsRB,focEhat);
                                        focEhat(amaxsave, &paramsRB);
                                        searchval_W = paramsRB.searchWoutST;
                                    }
                                }
                            }
                            valueEhat[inxRB(igridCOH,tgrid,zz)]  = -valfnmax;
                            C_Ehat[inxRB(igridCOH,tgrid,zz)]     = paramsRB.C_val;
                            saveEhat[inxRB(igridCOH,tgrid,zz)]   = amaxsave;
                            searchEhat[inxRB(igridCOH,tgrid,zz)] = searchval_W;
                            
                            
                            /** ----  Solving Bankruptcy case  ---- **/
                            icase = 100;
                            savemax=min(Gridmax-0.000000001,findconsRB((0.0),gridCOH[igridCOH])-0.000000001);
                            if(igridCOH == 0) {
                                valfnmax = focB(0.0,&paramsRB);
                                amaxsave = Gridmin; searchval_W = paramsRB.searchWoutST;
                                icase = 1.0;
                            } else {
                                if(focB(Gridmin,&paramsRB) < focB(Gridmin+0.0000001, &paramsRB)) {
                                    valfnmax = focB(0.0,&paramsRB);
                                    searchval_W = paramsRB.searchWoutST;  amaxsave = 0.0;
                                    icase = 1.0;
                                }else{//cas ou le max est a droite de la borne inf
                                    
                                    if(focB(savemax,&paramsRB) < focB((savemax-0.0000001), &paramsRB)) {
                                        valfnmax = focB(savemax,&paramsRB);
                                        searchval_W = paramsRB.searchWoutST; amaxsave = savemax;
                                        icase = 3.0;
                                    }
                                    
                                    if(icase >= 100) {
                                        valfnmax = mygolden(GridminCOH,(GridminCOH+0.0000001),savemax,TOL,amaxsave,&paramsRB,focB);
                                        focB(amaxsave, &paramsRB);
                                        searchval_W = paramsRB.searchWoutST;
                                    }
                                    
                                }
                            }
                            valueB[inxRB(igridCOH,tgrid,zz)]  = -valfnmax;
                            C_B[inxRB(igridCOH,tgrid,zz)]     = paramsRB.C_val;
                            saveB[inxRB(igridCOH,tgrid,zz)]   = amaxsave;
                            searchB[inxRB(igridCOH,tgrid,zz)] = searchval_W;
                            
                            
                            
                            /** ---- Solving REPAYMENT ---- **/
                            icase = 100;
                            savemax=min(Gridmax-0.000000001,findconsRB((0.0),gridCOH[igridCOH])-0.000000001);
                            if(igridCOH == 0) {
                                valfnmax = focR(0.0,&paramsRB);
                                searchval_W = paramsRB.searchWoutST; amaxsave = Gridmin;
                                icase = 1.0;
                            } else {
                                if(focR(0.0,&paramsRB) < focR(Gridmin+0.000001, &paramsRB)) {
                                    valfnmax = focR(0.0,&paramsRB);
                                    searchval_W = paramsRB.searchWoutST; amaxsave = Gridmin;
                                    icase = 1.0;
                                }else{//cas ou le max est a droite de la borne inf
                                    
                                    if(focR(savemax,&paramsRB) < focR((savemax-0.00001), &paramsRB)) {
                                        valfnmax = focR(savemax,&paramsRB);
                                        searchval_W = paramsRB.searchWoutST; amaxsave = savemax;
                                        icase = 3.0;
                                    }
                                    
                                    if(icase >= 100) {
                                        valfnmax = mygolden(Gridmin,(Gridmin+0.0000001),savemax,TOL,amaxsave,&paramsRB,focR);
                                        focR(amaxsave, &paramsRB);
                                        searchval_W = paramsRB.searchWoutST;
                                    }
                                    
                                }
                            }
                            valueR[inxRB(igridCOH, tgrid, zz)]  =  -valfnmax;
                            C_R[inxRB(igridCOH,tgrid,zz)]       =  paramsRB.C_val;
                            saveR[inxRB(igridCOH, tgrid, zz)]   =  amaxsave;
                            searchR[inxRB(igridCOH, tgrid, zz)] =  searchval_W;
                            
                            
                            // If Self-employment Assistance Program is activated.
#if SEA == 1
                            
                            paramsRB.expE_WWNE_JNEi_ST   = expE_WWNE_JNEi;
                            paramsRB.expE_WWE_JEi_ST     = expE_WWE_JEi;
                            paramsRB.expE_UiNE_JNEi_ST   = expE_UiNE_JNEi;
                            paramsRB.expE_UiE_JEi_ST     = expE_UiE_JEi;
                            paramsRB.expE_UiE_ST         = expE_UiE;
#if OPT_compute_C_flow == 1
                            paramsRB.C_expE_WWNE_JNEi_ST   = C_expE_WWNE_JNEi;
                            paramsRB.C_expE_WWE_JEi_ST     = C_expE_WWE_JEi;
                            paramsRB.C_expE_UiNE_JNEi_ST   = C_expE_UiNE_JNEi;
                            paramsRB.C_expE_UiE_JEi_ST     = C_expE_UiE_JEi;
                            paramsRB.C_expE_UiE_ST         = C_expE_UiE;
#endif
                            
                            paramsRB.JinsuredST          = 1;
                            
                            
                            /** ---- Solving Excluded HAT ---- **/
                            icase   =   100;
                            savemax =   min(Gridmax-0.000000001,findconsRB((0.0),gridCOH[igridCOH])-0.000000001);
                            if(igridCOH == 0) {
                                valfnmax = focEhat(0.0,&paramsRB); amaxsave = Gridmin; searchval_W = paramsRB.searchWoutST;
                                icase = 1.0;
                            } else {
                                if(focEhat(Gridmin,&paramsRB) < focEhat(Gridmin+0.0000001, &paramsRB)) {
                                    valfnmax = focEhat(0.0,&paramsRB);  searchval_W = paramsRB.searchWoutST;  amaxsave = 0.0;
                                    icase = 1.0;
                                }else{//cas ou le max est a droite de la borne inf
                                    
                                    if(focEhat(savemax,&paramsRB) < focEhat((savemax-0.0000001), &paramsRB)) {
                                        valfnmax = focEhat(savemax,&paramsRB); searchval_W = paramsRB.searchWoutST; amaxsave = savemax;
                                        icase = 3.0;
                                    }
                                    
                                    if(icase >= 100) {  // interior case.
                                        valfnmax = mygolden(GridminCOH,(GridminCOH+0.0000001),savemax,TOL,amaxsave,&paramsRB,focEhat);
                                        focEhat(amaxsave, &paramsRB);
                                        searchval_W = paramsRB.searchWoutST;
                                    }
                                }
                            }
                            valueEhati[inxRB(igridCOH,tgrid,zz)]  = -valfnmax;
                            C_Ehati[inxRB(igridCOH,tgrid,zz)]     = paramsRB.C_val;
                            saveEhati[inxRB(igridCOH,tgrid,zz)]   = amaxsave;
                            searchEhati[inxRB(igridCOH,tgrid,zz)] = searchval_W;
                            
                            
                            
                            /** ---- Solving BANKRUPT ---- **/
                            icase = 100;
                            savemax=min(Gridmax-0.000000001,findconsRB((0.0),gridCOH[igridCOH])-0.000000001);
                            if(igridCOH == 0) {
                                valfnmax = focB(0.0,&paramsRB);
                                amaxsave = Gridmin; searchval_W = paramsRB.searchWoutST;
                                icase = 1.0;
                            } else {
                                if(focB(Gridmin,&paramsRB) < focB(Gridmin+0.0000001, &paramsRB)) {
                                    valfnmax = focB(0.0,&paramsRB);
                                    searchval_W = paramsRB.searchWoutST;  amaxsave = 0.0;
                                    icase = 1.0;
                                }else{//cas ou le max est a droite de la borne inf
                                    
                                    if(focB(savemax,&paramsRB) < focB((savemax-0.0000001), &paramsRB)) {
                                        valfnmax = focB(savemax,&paramsRB);
                                        searchval_W = paramsRB.searchWoutST; amaxsave = savemax;
                                        icase = 3.0;
                                    }
                                    
                                    if(icase >= 100) {
                                        valfnmax = mygolden(GridminCOH,(GridminCOH+0.0000001),savemax,TOL,amaxsave,&paramsRB,focB);
                                        focB(amaxsave, &paramsRB);
                                        searchval_W = paramsRB.searchWoutST;
                                    }
                                    
                                }
                            }
                            valueBi[inxRB(igridCOH, tgrid, zz)]  = -valfnmax;
                            C_Bi[inxRB(igridCOH,tgrid,zz)]       = paramsRB.C_val;
                            saveBi[inxRB(igridCOH, tgrid, zz)]   = amaxsave;
                            searchBi[inxRB(igridCOH, tgrid, zz)] = searchval_W;
                            
                            
                            
                            /** ---- Solving REPAYMENT ---- **/
                            icase = 100;
                            savemax=min(Gridmax-0.000000001,findconsRB((0.0),gridCOH[igridCOH])-0.000000001);
                            if(igridCOH == 0) {
                                valfnmax = focR(0.0,&paramsRB);
                                searchval_W = paramsRB.searchWoutST; amaxsave = Gridmin;
                                icase = 1.0;
                            } else {
                                if(focR(0.0,&paramsRB) < focR(Gridmin+0.000001, &paramsRB)) {
                                    valfnmax = focR(0.0,&paramsRB);
                                    searchval_W = paramsRB.searchWoutST; amaxsave = Gridmin;
                                    icase = 1.0;
                                }else{//cas ou le max est a droite de la borne inf
                                    
                                    if(focR(savemax,&paramsRB) < focR((savemax-0.00001), &paramsRB)) {
                                        valfnmax = focR(savemax,&paramsRB);
                                        searchval_W = paramsRB.searchWoutST; amaxsave = savemax;
                                        icase = 3.0;
                                    }
                                    
                                    if(icase >= 100) {
                                        valfnmax = mygolden(Gridmin,(Gridmin+0.0000001),savemax,TOL,amaxsave,&paramsRB,focR);
                                        focR(amaxsave, &paramsRB);
                                        searchval_W = paramsRB.searchWoutST;
                                    }
                                    
                                }
                            }
                            valueRi[inxRB(igridCOH,tgrid,zz)]  =  -valfnmax;
                            C_Ri[inxRB(igridCOH,tgrid,zz)]     =  paramsRB.C_val;
                            saveRi[inxRB(igridCOH,tgrid,zz)]   =  amaxsave;
                            searchRi[inxRB(igridCOH,tgrid,zz)] =  searchval_W;
#endif
                            
                        } // end zz
                    } // end ygrid
                } // end igridCOH
#if OMP == 1
            } // end pragma
#endif
            
            
            
            
            /************************************************************************************/
            /** THIRD, solve the bankruptcy -- repayment, and rescale policy function.
             /************************************************************************************/
            
#if OMP == 1
#pragma omp parallel
            {
#pragma omp for reduction(max:critereVF) private(params,tgrid, igrid, zgrid, zz, bb, znext, ygrid, igridCOH, ygridnext, tgridnext, zgridnext, iexcluded, ilongrun, iinsured, icaseJ, icase, ixgridB, ixgridR, valmax, bruteax, amaxAc, savemax, amaxsave, searchval_J, searchval_W, valfnmax, cohB, cohR, xgridB, dxgridB, xgridR, dxgridR, Bval, Rval, focE_grid, focE_gridm01, focE_min, focE_minp01, focE_max, focE_maxm01, RB_max,valopt1, axopti1, rate_endo_opti,valopt2, axopti2, lowerbound, upperbound,valopt3, axopti3, pointbb3_aft, pointbb3_bef,valueJNEtemp, saveJNEtemp, searchJNEtemp, defaultJNEtemp, valueJEtemp, saveJEtemp, searchJEtemp,tempE, COHB, COHR, valueBB, valueRR, profitR, muENDO,RATEendo,collateral,C_tempE,C_BB,C_RR)
#endif
                for(igrid=0;igrid<maxgrid;igrid++){
                    for(tgrid=0;tgrid<maxtheta;tgrid++){
                        
                        params.iterST=iter;
                        params.igridST=igrid; params.tgridST=tgrid;
                        
                        params.valueBnextST     =   valueB;
                        params.valueRnextST     =   valueR;
#if OPT_compute_C_flow == 1
                        params.C_BnextST        =   C_B;
                        params.C_RnextST        =   C_R;
                        params.C_EhatnextST     =   C_Ehat;
#endif
                        params.saveBnextST      =   saveB;
                        params.saveRnextST      =   saveR;
                        params.searchBnextST    =   searchB;
                        params.searchRnextST    =   searchR;
                        params.valueEhatnextST  =   valueEhat;
                        params.saveEhatnextST   =   saveEhat;
                        params.searchEhatnextST =   searchEhat;
                        
#if SEA == 1
                        params.valueBinextST    =   valueBi;
                        params.valueRinextST    =   valueRi;
#if OPT_compute_C_flow == 1
                        params.C_BinextST       =   C_Bi;
                        params.C_RinextST       =   C_Ri;
                        params.C_EhatinextST    =   C_Ehati;
#endif
                        params.saveBinextST     =   saveBi;
                        params.saveRinextST     =   saveRi;
                        params.searchBinextST   =   searchBi;
                        params.searchRinextST   =   searchRi;
                        params.valueEhatinextST =   valueEhati;
                        params.saveEhatinextST  =   saveEhati;
                        params.searchEhatinextST=   searchEhati;
#endif
                        
                        /** ----- Solving the ENTREPRENEUR problem ----- **/
                        // -- solve the capital problem --
                        
                        for(zgrid=0;zgrid<maxfirmtype;zgrid++) {
                            
                            
#if SEA == 1
                            for(iinsured = 0; iinsured < 2; iinsured++){
#endif
                                
#if SEA == 0
                                iinsured = 0;
#endif
                                
                                params.JinsuredST   = iinsured;
                                params.zgridST      = zgrid;
                                
                                
                                /***************************/
                                /** ENTREPRENEUR  NOT EXCLUDED **/
                                /***************************/
                                
                                
                                // apply the howard improvement :: compute the optimal rate + capital investment only every X iteration.
                                
                                if(iter < 5 || howard_off == 1.0 || iter % n_iter == 0){
                                    
                                    // bracketing ---
                                    focE_max        = focE(Gridmax+1000,&params);
                                    focE_maxm01     = focE((Gridmax+1000-0.0000001),&params);
                                    focE_min        = focE(0.0,&params);
                                    //focE_minp01     = focE(0.0000001,&params);
                                    
                                    // Strategy to compute the entrepreneur's problem
                                    // Golden between 0 to ax_threshold, then grid search
                                    icase = 100;
                                    
                                    if(igrid == 0) { // cannot invest sufficiently
                                        valfnmax = focE_min;    amaxAc = 0.0;   icase = 0.0;
                                    }else{
                                        //if(focE_min < focE_minp01) { // check corner solution at the bottom
                                        //    valfnmax = focE_min;    amaxAc = 0.0;   icase = 1.0;
                                        //} else {

                                            // multi-search minimum //
                                            valopt1 = 1000000; valopt2 = 1000000; valopt3 = 1000000;
                                            
                                            /** METHOD 1, between [0, a] :: CAN BE REMOVED LATER **/
                                            params.methodST = 1;
                                            focE_grid       = focE((grid[igrid]),&params);
                                            focE_gridm01    = focE(grid[igrid]-0.0000001,&params);
                                            
                                            if(focE_gridm01 < focE_grid) {
                                                valopt1 = mygolden(0.0,0.0000001,grid[igrid],TOL,axopti1,&params,focE);
//                                                icase   = 110;
                                            } else {
                                                valopt1 = focE_grid;
                                                axopti1 = grid[igrid];
                                            }
                                            
//                                            if(icase == 110){  // end icase == 110 where minimum is before between [0, a]
//                                                amaxAc = axopti1;   valfnmax = valopt1;
//                                            }else{
                                                
                                                /** METHOD 2, between [0, axmax]  -- simple goldensearch **/
                                                params.methodST = 2;
                                                if(focE_max < focE_maxm01) {
                                                    valopt2 = focE_max;
                                                    axopti2 = Gridmax+1000;
                                                }else{
                                                    valopt2 = mygolden(0.00000001,grid[igrid],Gridmax+1000,TOL,axopti2,&params,focE);
                                                }
                                                
                                                /** METHOD 3, between [axz0, (1+lambda)*a] apply a mix gridsearch - golden **/
                                                params.methodST = 3;
                                                lowerbound = 0.0; //grid[igrid];
                                                upperbound = Gridmax+1000;//;grid[igrid]*(1.0+lambdapar);
                                                
                                                // I:: BRACKET the minimum using a gridsearch //
                                                for(bb=0;bb<maxAX;bb++){
                                                    bruteax  = expspace(bb,lowerbound,upperbound,Echelle1,(maxAX));
                                                    valfnmax = focE(bruteax, &params);
                                                    if(valfnmax < valopt3){
                                                        valopt3 = focE(bruteax, &params);
                                                        axopti3 = bruteax;
                                                        if(bb > (0)){
                                                            tmp_point = expspace((bb-1),lowerbound,upperbound,Echelle1,(maxAX));
                                                            pointbb3_bef = max(tmp_point,lowerbound);
                                                        }else{
                                                            pointbb3_bef = lowerbound;
                                                        }
                                                        if(bb < (maxAX-1)){
                                                            pointbb3_aft = min(expspace((bb+1),lowerbound,upperbound,Echelle1,(maxAX)),upperbound);
                                                        }else{
                                                            pointbb3_aft = upperbound;
                                                        }
                                                    }
                                                }
                                        
                                                // II:: GOLDEN and find the minimum using previous bracketing //
                                                valopt3 = mygolden(pointbb3_bef,pointbb3_bef+0.0000001,pointbb3_aft,TOL,axopti3,&params,focE);
                                                valopt3 = focE(axopti3, &params);

                                                /** LAST STEP:: choose the best solution **/
                                                if(valopt2 <= valopt1 && valopt2 <= valopt3){amaxAc = axopti2;  valfnmax = valopt2;}
                                                if(valopt1 <= valopt2 && valopt1 <= valopt3){amaxAc = axopti1;  valfnmax = valopt1;} // printf("ERROR? %f %f %f | %f %f %f ",valopt1,valopt2,valopt3,axopti1,axopti2,axopti3);getchar();}
                                                if(valopt3 <= valopt2 && valopt3 <= valopt1){amaxAc = axopti3;  valfnmax = valopt3;}
                                                
                                                //amaxAc = axopti3;  valfnmax = valopt3;
                                                // printf("checks :: %f %f %f", lowerbound,upperbound,axopti3);getchar();
                                                
//                                            }
                                            
//                                        } // end ax > 0.0000000001
                                    } // end igrid == 0 (case where you can not invest)
                                    
                                    
                                    focE(amaxAc,&params); // to make sure we have the right params.
                                    rate_endo_opti = params.rate_endo;
                                    
                                    
                                    /** UNINSURED ENTREPRENEUR **/
                                    if(iinsured == 0){
                                        critereVF = max(critereVF,fabs(valueJNE[inxE(igrid, tgrid, zgrid)] - (-valfnmax)));
                                        collateralJNE[inxE(igrid, tgrid, zgrid)]    =   amaxAc;
                                        if(amaxAc>grid[igrid]){sloanJNE[inxE(igrid,tgrid,zgrid)]=(amaxAc-grid[igrid]);}else{sloanJNE[inxE(igrid,tgrid,zgrid)]=0.0;}
                                        valueJNE[inxE(igrid, tgrid, zgrid)]         =   -valfnmax;
                                        C_JNE[inxE(igrid, tgrid, zgrid)]            =   params.C_val;
                                        rateJNE[inxE(igrid, tgrid, zgrid)]          =   rate_endo_opti;
                                        
                                        if(igrid == 0) {
                                            for(zz=0; zz<maxfirmtype; zz++) {
                                                valueJNEz[inxEE(igrid, tgrid, zgrid, zz)]   = -100000.0;
                                                saveJNE[inxEE(igrid, tgrid, zgrid, zz)]     = 0.0;
                                                defaultJNE[inxEE(igrid, tgrid, zgrid, zz)]  = 1.0;
                                                searchJNE_W[inxEE(igrid, tgrid, zgrid, zz)] = 0.0;
                                            }
                                        }else{
                                            
                                            valueJNEtemp    = params.JNE_val;
                                            saveJNEtemp     = params.JNE_save;
                                            searchJNEtemp   = params.JNE_search;
                                            defaultJNEtemp  = params.JNE_default;
                                            
                                            for(zz=0; zz<maxfirmtype; zz++) {
                                                valueJNEz[inxEE(igrid, tgrid, zgrid, zz)]   = valueJNEtemp[zz];
                                                saveJNE[inxEE(igrid, tgrid, zgrid, zz)]     = saveJNEtemp[zz];
                                                defaultJNE[inxEE(igrid, tgrid, zgrid, zz)]  = defaultJNEtemp[zz];
                                                searchJNE_W[inxEE(igrid, tgrid, zgrid, zz)] = searchJNEtemp[zz];
                                            }
                                        }
                                    }
                                    
                                    /** INSURED ENTREPRENEUR **/
#if SEA == 1
                                    if(iinsured == 1){
                                        critereVF = max(critereVF,fabs(valueJNEi[inxE(igrid, tgrid, zgrid)] - (-valfnmax)));
                                        collateralJNEi[inxE(igrid, tgrid, zgrid)]   =   amaxAc;
                                        if(amaxAc>grid[igrid]){sloanJNEi[inxE(igrid, tgrid, zgrid)]=(amaxAc-grid[igrid]);}else{sloanJNEi[inxE(igrid, tgrid, zgrid)]=0.0;}
                                        valueJNEi[inxE(igrid, tgrid, zgrid)]        =   -valfnmax;
                                        C_JNEi[inxE(igrid, tgrid, zgrid)]           =   params.C_val;
                                        rateJNEi[inxE(igrid, tgrid, zgrid)]         =   rate_endo_opti;
                                        
                                        if(igrid == 0) {
                                            for(zz=0; zz<maxfirmtype; zz++) {
                                                valueJNEiz[inxEE(igrid, tgrid, zgrid, zz)]      = -100000.0;
                                                saveJNEi[inxEE(igrid, tgrid, zgrid, zz)]        = 0.0;
                                                defaultJNEi[inxEE(igrid, tgrid, zgrid, zz)]     = 1.0;
                                                searchJNEi_W[inxEE(igrid, tgrid, zgrid, zz)]    = 0.0;
                                            }
                                        } else {
                                            valueJNEtemp    = params.JNE_val;
                                            saveJNEtemp     = params.JNE_save;
                                            searchJNEtemp   = params.JNE_search;
                                            defaultJNEtemp  = params.JNE_default;
                                            
                                            for(zz=0; zz<maxfirmtype; zz++) {
                                                valueJNEiz[inxEE(igrid, tgrid, zgrid, zz)]      = valueJNEtemp[zz];
                                                saveJNEi[inxEE(igrid, tgrid, zgrid, zz)]        = saveJNEtemp[zz];
                                                defaultJNEi[inxEE(igrid, tgrid, zgrid, zz)]     = defaultJNEtemp[zz];
                                                searchJNEi_W[inxEE(igrid, tgrid, zgrid, zz)]    = searchJNEtemp[zz];
                                            }
                                        }
                                    }
#endif
                                    
                                }else{
                                    
                                    // bank condition first guess //
                                    tempE = 0.0; C_tempE = 0.0;
                                    if(iinsured == 0){
                                        RATEendo   = rateJNE[inxE(igrid, tgrid, zgrid)];
                                        collateral = collateralJNE[inxE(igrid, tgrid, zgrid)];
                                    }
#if SEA == 1
                                    if(iinsured == 1){
                                        RATEendo   = rateJNEi[inxE(igrid, tgrid, zgrid)];
                                        collateral = collateralJNEi[inxE(igrid, tgrid, zgrid)];
                                    }
#endif
                                    
                                    for(zz = 0; zz<maxfirmtype; zz++){
                                        
                                        COHB    = findCOH_B(igrid, zz, tgrid, collateral, RATEendo, iinsured);
                                        xgridB  = invexpspace(COHB,GridminCOH,GridmaxCOH,EchelleCOH,maxgridCOH);
                                        dxgridB = weightinter(COHB, xgridB, &ixgridB, gridCOH, maxgridCOH);
                                        
                                        COHR    = max(0.0,findCOH_R(igrid, zz, tgrid, collateral, RATEendo, iinsured));
                                        xgridR  = invexpspace(COHR,GridminCOH,GridmaxCOH,EchelleCOH,maxgridCOH);
                                        dxgridR = weightinter(COHR, xgridR, &ixgridR, gridCOH, maxgridCOH);
                                        
                                        // COMPUTE THE INDUCED POLICIES -- normal case.
                                        if(iinsured == 0){
                                            valueBB  = inter1d(dxgridB,valueB[inxRB((ixgridB),tgrid,zz)],valueB[inxRB((ixgridB+1),tgrid,zz)]);
                                            valueRR  = inter1d(dxgridR,valueR[inxRB((ixgridR),tgrid,zz)],valueR[inxRB((ixgridR+1),tgrid,zz)]);
                                            C_BB  = inter1d(dxgridB,C_B[inxRB((ixgridB),tgrid,zz)],C_B[inxRB((ixgridB+1),tgrid,zz)]);
                                            C_RR  = inter1d(dxgridR,C_R[inxRB((ixgridR),tgrid,zz)],C_R[inxRB((ixgridR+1),tgrid,zz)]);
                                        }
#if SEA == 1
                                        if(iinsured == 1){
                                            muENDO  = muendoSEA(0.0,tgrid); //mupar*(bstarcostfun(0.0,tgrid)/((1.0-taxrate)*rhostar*wstar*prod[tgrid]));
                                            
                                            valueBB  = muENDO*inter1d(dxgridB,valueB[inxRB((ixgridB),tgrid,zz)],valueB[inxRB((ixgridB+1),tgrid,zz)]) + (1.0-muENDO)*inter1d(dxgridB,valueBi[inxRB((ixgridB),tgrid,zz)],valueBi[inxRB((ixgridB+1),tgrid,zz)]);
                                            C_BB  = muENDO*inter1d(dxgridB,C_B[inxRB((ixgridB),tgrid,zz)],C_B[inxRB((ixgridB+1),tgrid,zz)]) + (1.0-muENDO)*inter1d(dxgridB,C_Bi[inxRB((ixgridB),tgrid,zz)],C_Bi[inxRB((ixgridB+1),tgrid,zz)]);
                                            
                                            profitR = profitRfun(igrid, zz, tgrid, collateral, RATEendo);
                                            muENDO  = muendoSEA(profitR,tgrid); //mupar*(bstarcostfun(profitR,tgrid)/((1.0-taxrate)*rhostar*wstar*prod[tgrid]));
                                            
                                            valueRR  = muENDO*inter1d(dxgridR,valueR[inxRB((ixgridR),tgrid,zz)],valueR[inxRB((ixgridR+1),tgrid,zz)]) + (1.0-muENDO)*inter1d(dxgridR,valueRi[inxRB((ixgridR),tgrid,zz)],valueRi[inxRB((ixgridR+1),tgrid,zz)]);
                                            C_RR  = muENDO*inter1d(dxgridR,C_R[inxRB((ixgridR),tgrid,zz)],C_R[inxRB((ixgridR+1),tgrid,zz)]) + (1.0-muENDO)*inter1d(dxgridR,C_Ri[inxRB((ixgridR),tgrid,zz)],C_Ri[inxRB((ixgridR+1),tgrid,zz)]);
                                        }
#endif
                                        
                                        if((valueBB > valueRR) || findCOH_R(igrid, zz, tgrid, collateral, RATEendo, iinsured) < 0.0){
                                            tempE   += mzprob[zgrid][zz]*valueBB;
                                            C_tempE += mzprob[zgrid][zz]*C_BB;
                                        }
                                        if((valueBB <= valueRR)){
                                            tempE   += mzprob[zgrid][zz]*valueRR;
                                            C_tempE += mzprob[zgrid][zz]*C_RR;
                                        }
                                    }
                                    
                                    if(iinsured == 0){
                                        critereVF = max(critereVF,fabs(valueJNE[inxE(igrid, tgrid, zgrid)] - tempE));
                                        valueJNE[inxE(igrid, tgrid, zgrid)] = tempE;
                                        C_JNE[inxE(igrid, tgrid, zgrid)]    = C_tempE;
                                    }
#if SEA == 1
                                    if(iinsured == 1){
                                        critereVF = max(critereVF,fabs(valueJNEi[inxE(igrid, tgrid, zgrid)] - tempE));
                                        valueJNEi[inxE(igrid, tgrid, zgrid)] = tempE;
                                        C_JNEi[inxE(igrid, tgrid, zgrid)]    = C_tempE;
                                    }
#endif
                                    
                                }
                                
                                
                                
                                /***************************/
                                /** ENTREPRENEUR EXCLUDED **/
                                /***************************/
                                if(iter < 5 || howard_off == 1.0  || iter % 2 == 0){
                                    
//                                    icaseJ = 100;
//                                    if(igrid == 0) { // cannot invest
//                                        valfnmax = focEC(0.0,&params);  amaxAc = 0.0;   icaseJ = 0.0;
//                                    }else{
//                                        if(focEC(0.0,&params) < focEC(TOL, &params)){valfnmax = focEC(0.0,&params); amaxAc = 0.0; icaseJ = 1.0;}
//                                        if(focEC(0.0,&params) >= focEC(TOL, &params)) {
//                                            if(focEC(grid[igrid],&params) < focEC(grid[igrid]-TOL, &params)) { // corner solution at the top
//                                                valfnmax = focEC(grid[igrid],&params); amaxAc = grid[igrid]; icaseJ = 3.0;
//                                            }
//                                        }
//                                        if(icaseJ >= 100) { // interior solution
//                                            valfnmax = mygolden(Gridmin,(Gridmin+TOL),grid[igrid],TOL,amaxAc,&params,focEC);
//                                            focEC(amaxAc, &params);
//                                        }
//                                    }
                                    
                                    // bracketing ---
                                    focE_max        = focEC(grid[igrid],&params);
                                    focE_maxm01     = focEC((grid[igrid]-0.0000001),&params);
                                    focE_min        = focEC(0.0,&params);
                                    //focE_minp01     = focE(0.0000001,&params);
                                    
                                    // Strategy to compute the entrepreneur's problem
                                    // Golden between 0 to ax_threshold, then grid search
                                    icase = 100;
                                    
                                    if(igrid == 0) { // cannot invest sufficiently
                                        valfnmax = focEC(0.0,&params);  amaxAc = 0.0;   icaseJ = 0.0;
                                    }else{
                                        //if(focE_min < focE_minp01) { // check corner solution at the bottom
                                        //    valfnmax = focE_min;    amaxAc = 0.0;   icase = 1.0;
                                        //} else {
                                        
                                        // multi-search minimum //
                                        valopt1 = 1000000; valopt2 = 1000000; valopt3 = 1000000;
                                        
                                        /** METHOD 1, between [0, a] :: CAN BE REMOVED LATER **/
                                        params.methodST = 1;
                                        focE_grid       = focEC((grid[igrid]),&params);
                                        focE_gridm01    = focEC(grid[igrid]-0.0000001,&params);
                                        
                                        if(focE_gridm01 < focE_grid) {
                                            valopt1 = mygolden(0.0,grid[igrid]/2,grid[igrid],TOL,axopti1,&params,focEC);
                                            //                                                icase   = 110;
                                        } else {
                                            valopt1 = focE_grid;
                                            axopti1 = grid[igrid];
                                        }
                                        
                                        //                                            if(icase == 110){  // end icase == 110 where minimum is before between [0, a]
                                        //                                                amaxAc = axopti1;   valfnmax = valopt1;
                                        //                                            }else{
                                        
//                                        /** METHOD 2, between [0, axmax]  -- simple goldensearch **/
//                                        params.methodST = 2;
//                                        if(focE_max < focE_maxm01) {
//                                            valopt2 = focE_max;
//                                            axopti2 = grid[igrid];
//                                        }else{
//                                            valopt2 = mygolden(0.00000001,grid[igrid]/2,grid[igrid],TOL,axopti2,&params,focEC);
//                                        }
                                        
                                        /** METHOD 3, between [axz0, (1+lambda)*a] apply a mix gridsearch - golden **/
                                        params.methodST = 3;
                                        lowerbound = 0.0; //grid[igrid];
                                        upperbound = grid[igrid];//;grid[igrid]*(1.0+lambdapar);
                                        
                                        // I:: BRACKET the minimum using a gridsearch //
                                        for(bb=0;bb<maxAX;bb++){
                                            bruteax  = expspace(bb,lowerbound,upperbound,Echelle1,(maxAX));
                                            valfnmax = focEC(bruteax, &params);
                                            if(valfnmax < valopt3){
                                                valopt3 = focEC(bruteax, &params);
                                                axopti3 = bruteax;
                                                if(bb > (0)){
                                                    tmp_point = expspace((bb-1),lowerbound,upperbound,Echelle1,(maxAX));
                                                    pointbb3_bef = max(tmp_point,lowerbound);
                                                }else{
                                                    pointbb3_bef = lowerbound;
                                                }
                                                if(bb < (maxAX-1)){
                                                    pointbb3_aft = min(expspace((bb+1),lowerbound,upperbound,Echelle1,(maxAX)),upperbound);
                                                }else{
                                                    pointbb3_aft = upperbound;
                                                }
                                            }
                                        }
                                        
                                        // II:: GOLDEN and find the minimum using previous bracketing //
                                        valopt3 = mygolden(pointbb3_bef,pointbb3_bef+0.0000001,pointbb3_aft,TOL,axopti3,&params,focEC);
                                        valopt3 = focEC(axopti3, &params);
                                        
                                        /** LAST STEP:: choose the best solution **/
                                        //if(valopt2 <= valopt1 && valopt2 <= valopt3){amaxAc = axopti2;  valfnmax = valopt2;}
                                        if(valopt1 <= valopt2 && valopt1 <= valopt3){amaxAc = axopti1;  valfnmax = valopt1;} // printf("ERROR? %f %f %f | %f %f %f ",valopt1,valopt2,valopt3,axopti1,axopti2,axopti3);getchar();}
                                        if(valopt3 <= valopt2 && valopt3 <= valopt1){amaxAc = axopti3;  valfnmax = valopt3;}
                                        
                                        //amaxAc = axopti3;  valfnmax = valopt3;
                                        // printf("checks :: %f %f %f", lowerbound,upperbound,axopti3);getchar();
                                        
                                        //                                            }
                                        
                                        //                                        } // end ax > 0.0000000001
                                    }
                                }else{
                                    
                                    if(iinsured == 0){
                                        valfnmax = focEC(collateralJE[inxE(igrid, tgrid, zgrid)],&params);
                                        amaxAc   = collateralJE[inxE(igrid, tgrid, zgrid)];
                                    }
                                    
#if SEA == 1
                                    if(iinsured == 1){
                                        valfnmax = focEC(collateralJEi[inxE(igrid, tgrid, zgrid)],&params);
                                        amaxAc   = collateralJEi[inxE(igrid, tgrid, zgrid)];
                                    }
#endif
                                }
                                
                                /** NOT INSURED ENTREPRENEUR **/
                                if(iinsured == 0){
                                    if(iter < 5 || iter % n_iter == 0){critereVF = max(critereVF,fabs(valueJE[inxE(igrid, tgrid, zgrid)] - (-valfnmax)));}
                                    collateralJE[inxE(igrid, tgrid, zgrid)] =    amaxAc;
                                    valueJE[inxE(igrid, tgrid, zgrid)]      =   -valfnmax;
                                    C_JE[inxE(igrid, tgrid, zgrid)]         =   params.C_val;
                                    
                                    if(igrid == 0) {
                                        for(zz=0; zz<maxfirmtype; zz++) {
                                            valueJEz[inxEE(igrid, tgrid, zgrid, zz)]    = -100000.0;
                                            saveJE[inxEE(igrid, tgrid, zgrid, zz)]      = 0.0;
                                            searchJE_W[inxEE(igrid, tgrid, zgrid, zz)]  = 0.0;
                                        }
                                    } else {
                                        valueJEtemp     = params.JE_val;
                                        saveJEtemp      = params.JE_save;
                                        searchJEtemp    = params.JE_search;
                                        
                                        for(zz=0; zz<maxfirmtype; zz++) {
                                            valueJEz[inxEE(igrid, tgrid, zgrid, zz)]    = valueJEtemp[zz];
                                            saveJE[inxEE(igrid, tgrid, zgrid, zz)]      = saveJEtemp[zz];
                                            searchJE_W[inxEE(igrid, tgrid, zgrid, zz)]  = searchJEtemp[zz];
                                        }
                                    }
                                }
                                
                                /** INSURED ENTREPRENEUR **/
#if SEA == 1
                                if(iinsured == 1){
                                    if(iter < 5 || iter % n_iter == 0){critereVF = max(critereVF,fabs(valueJEi[inxE(igrid, tgrid, zgrid)] - (-valfnmax)));}
                                    collateralJEi[inxE(igrid, tgrid, zgrid)]    =   amaxAc;
                                    valueJEi[inxE(igrid, tgrid, zgrid)]         =   -valfnmax;
                                    C_JEi[inxE(igrid, tgrid, zgrid)]            =   params.C_val;
                                    
                                    if(igrid == 0) {
                                        for(zz=0; zz<maxfirmtype; zz++) {
                                            valueJEiz[inxEE(igrid, tgrid, zgrid, zz)]   = -100000.0;
                                            saveJEi[inxEE(igrid, tgrid, zgrid, zz)]     = 0.0;
                                            searchJEi_W[inxEE(igrid, tgrid, zgrid, zz)] = 0.0;
                                        }
                                    } else {
                                        valueJEtemp     = params.JE_val;
                                        saveJEtemp      = params.JE_save;
                                        searchJEtemp    = params.JE_search;
                                        
                                        for(zz=0; zz<maxfirmtype; zz++) {
                                            valueJEiz[inxEE(igrid, tgrid, zgrid, zz)]   = valueJEtemp[zz];
                                            saveJEi[inxEE(igrid, tgrid, zgrid, zz)]     = saveJEtemp[zz];
                                            searchJEi_W[inxEE(igrid, tgrid, zgrid, zz)] = searchJEtemp[zz];
                                        }
                                    }
                                }
#endif
                                
#if SEA == 1
                            } // loop over insurance state variable
#endif
                            
                            
                            
                        } // end zgrid
                        
                        
                        
                        
                        
                        /***********************************/
                        /** Solving the UNEMPLOYED problem //
                         /***********************************/
                        
                        params.expU_UnNE_JNE_ST      = expU_UnNE_JNE;
                        params.expU_UiNE_JNE_ST      = expU_UiNE_JNE;
                        params.expU_UnE_JE_ST        = expU_UnE_JE;
                        params.expU_UiE_JE_ST        = expU_UiE_JE;
                        
                        params.expU_UnNE_ST          = expU_UnNE;
                        params.expU_UnE_ST           = expU_UnE;
                        params.expU_UiNE_ST          = expU_UiNE;
                        params.expU_UiE_ST           = expU_UiE;
                        
                        params.expU_WWNE_ST          = expU_WWNE;
                        params.expU_WWE_ST           = expU_WWE;
                        params.expU_WWNE_JNE_ST      = expU_WWNE_JNE;
                        params.expU_WWE_JE_ST        = expU_WWE_JE;
#if OPT_compute_C_flow == 1
                        params.C_expU_UnNE_JNE_ST      = C_expU_UnNE_JNE;
                        params.C_expU_UiNE_JNE_ST      = C_expU_UiNE_JNE;
                        params.C_expU_UnE_JE_ST        = C_expU_UnE_JE;
                        params.C_expU_UiE_JE_ST        = C_expU_UiE_JE;
                        
                        params.C_expU_UnNE_ST          = C_expU_UnNE;
                        params.C_expU_UnE_ST           = C_expU_UnE;
                        params.C_expU_UiNE_ST          = C_expU_UiNE;
                        params.C_expU_UiE_ST           = C_expU_UiE;
                        
                        params.C_expU_WWNE_ST          = C_expU_WWNE;
                        params.C_expU_WWE_ST           = C_expU_WWE;
                        params.C_expU_WWNE_JNE_ST      = C_expU_WWNE_JNE;
                        params.C_expU_WWE_JE_ST        = C_expU_WWE_JE;
#endif
                        
#if SEA == 1
                        params.expU_UiNE_JNEi_ST     = expU_UiNE_JNEi;
                        params.expU_UiE_JEi_ST       = expU_UiE_JEi;
                        params.expU_WWNE_JNEi_ST     = expU_WWNE_JNEi;
                        params.expU_WWE_JEi_ST       = expU_WWE_JEi;
#if OPT_compute_C_flow == 1
                        params.C_expU_UiNE_JNEi_ST     = C_expU_UiNE_JNEi;
                        params.C_expU_UiE_JEi_ST       = C_expU_UiE_JEi;
                        params.C_expU_WWNE_JNEi_ST     = C_expU_WWNE_JNEi;
                        params.C_expU_WWE_JEi_ST       = C_expU_WWE_JEi;
#endif
#endif
                        
                        for(iexcluded = 0; iexcluded < 2; iexcluded++){
                            for(ilongrun = 0; ilongrun < 2; ilongrun++){
                                
                                params.excludedST   = iexcluded;    // iexcluded == 1 ** IF EXCLUDED **   iexcluded == 0 ** IF NOT EXCLUDED **
                                params.Utype        = ilongrun;     // ilongrun == 0 ** IF SHORTRUN U **  ilongrun == 1  ** IF LONGRUN U **
                                
                                if (focU_E(Gridmin,&params)<focU_E(Gridmin+TOL,&params)){// corner solution at the bottom
                                    valfnmax    =   focU_E(Gridmin,&params);
                                    amaxsave    =   Gridmin;
                                    searchval_W =   params.searchWoutST; searchval_J =   params.searchJoutST;
                                }else{//cas ou le max est a droite de la borne inf
                                    
                                    savemax = min(Gridmax-0.000000001,findconsU((0.0),igrid, tgrid, ilongrun)-0.000000001);
                                    
                                    // Si utility(today avec un peu moins) < utility(today avec un peu plus) alors intérêt encore à épargner (monotonicity).
                                    // Recall that focWW give the opposite of the valuefunction (find min)
                                    if (focU_E(savemax-TOL,&params)>focU_E(savemax,&params)){ // corner solution at the top.
                                        valfnmax    =   focU_E(savemax, &params);
                                        searchval_W =   params.searchWoutST; searchval_J = params.searchJoutST;
                                        amaxsave    =   savemax;
                                    }else{// interior solution
                                        valfnmax    =   mygolden(Gridmin,(Gridmin+TOL),savemax,TOL,amaxsave,&params,focU_E);
                                        valfnmax    =   focU_E(amaxsave, &params);
                                        searchval_W =   params.searchWoutST; searchval_J = params.searchJoutST;
                                    }
                                }
                                
                                if(iexcluded == 1 && ilongrun == 0){
                                    critereVF = max(critereVF,fabs(valueUiE[inx(igrid, tgrid)] - (-valfnmax)));
                                    searchUiE_W[inx(igrid, tgrid)] = searchval_W;   searchUiE_J[inx(igrid, tgrid)] = searchval_J;
                                    valueUiE[inx(igrid, tgrid)] = -valfnmax; 
                                    C_UiE[inx(igrid, tgrid)]    = params.C_val;
                                    saveUiE[inx(igrid, tgrid)]  = amaxsave;
                                }
                                
                                if(iexcluded == 0 && ilongrun == 0){
                                    critereVF = max(critereVF,fabs(valueUiNE[inx(igrid, tgrid)] - (-valfnmax)));
                                    searchUiNE_W[inx(igrid, tgrid)] = searchval_W;  searchUiNE_J[inx(igrid, tgrid)] = searchval_J;
                                    valueUiNE[inx(igrid, tgrid)] = -valfnmax;
                                    C_UiNE[inx(igrid, tgrid)]    = params.C_val;
                                    saveUiNE[inx(igrid, tgrid)] = amaxsave;
                                }
                                
                                if(iexcluded == 1 && ilongrun == 1){
                                    critereVF = max(critereVF,fabs(valueUnE[inx(igrid, tgrid)] - (-valfnmax)));
                                    searchUnE_W[inx(igrid, tgrid)] = searchval_W;   searchUnE_J[inx(igrid, tgrid)] = searchval_J;
                                    valueUnE[inx(igrid, tgrid)] = -valfnmax;
                                    C_UnE[inx(igrid, tgrid)]    = params.C_val;
                                    saveUnE[inx(igrid, tgrid)]  = amaxsave;
                                }
                                
                                if(iexcluded == 0 && ilongrun == 1){
                                    critereVF = max(critereVF,fabs(valueUnNE[inx(igrid, tgrid)] - (-valfnmax)));
                                    searchUnNE_W[inx(igrid, tgrid)] = searchval_W;  searchUnNE_J[inx(igrid, tgrid)] = searchval_J;
                                    valueUnNE[inx(igrid, tgrid)] = -valfnmax;    
                                    C_UnNE[inx(igrid, tgrid)]    = params.C_val;
                                    saveUnNE[inx(igrid, tgrid)]  = amaxsave;
                                }
                                
                            } // end ishortrun
                        } // end iexcluded
                        
                        
                        
                        
                        
                        /********************************/
                        /** Solving the WORKER problem **/
                        /********************************/
                        
                        for(ygrid = 0; ygrid < maxyprod; ygrid++) {     // boucle for ygrid values.
                            
                            params.ygridST = ygrid;
                            
                            params.expW_WWE_JE_ST          = expW_WWE_JE;
                            params.expW_WWNE_JNE_ST        = expW_WWNE_JNE;
                            params.expW_WWE_ST             = expW_WWE;
                            params.expW_WWNE_ST            = expW_WWNE;
#if OPT_compute_C_flow == 1
                            params.C_expW_WWE_JE_ST          = C_expW_WWE_JE;
                            params.C_expW_WWNE_JNE_ST        = C_expW_WWNE_JNE;
                            params.C_expW_WWE_ST             = C_expW_WWE;
                            params.C_expW_WWNE_ST            = C_expW_WWNE;
#endif
                            
                            for(iexcluded = 0; iexcluded < 2; iexcluded++){
                                
                                params.excludedST = iexcluded;
                                
                                if(focW_E(Gridmin,&params)<focW_E(Gridmin+TOL,&params)){// corner solution at the bottom
                                    valfnmax = focW_E(Gridmin,&params);     amaxsave=Gridmin;   searchval_J = params.searchJoutST;
                                }else{//cas ou le max est a droite de la borne inf
                                    
                                    // define the highest level of xval, which means that consume 0.000000001 and kept the rest for saving.
                                    savemax = min(Gridmax-0.000000001,findconsW((0.0),igrid, tgrid, ygrid)-0.000000001);
                                    
                                    if(focW_E(savemax-TOL,&params)>focW_E(savemax,&params)){ // corner solution at the top.
                                        valfnmax    =   focW_E(savemax, &params);
                                        amaxsave    =   savemax;
                                        searchval_J =   params.searchJoutST;
                                    }else{// interior solution
                                        valfnmax    =   mygolden(Gridmin,(Gridmin+TOL),savemax,TOL,amaxsave,&params,focW_E);
                                        focW_E(amaxsave, &params);
                                        searchval_J =   params.searchJoutST;
                                    }
                                }
                                
                                if(iexcluded == 1){
                                    critereVF = max(critereVF,fabs(valueWWE[inxW(igrid, tgrid, ygrid)] - (-valfnmax)));
                                    valueWWE[inxW(igrid, tgrid, ygrid)]     =   -valfnmax;
                                    C_WWE[inxW(igrid, tgrid, ygrid)]        =   params.C_val;
                                    saveWWE[inxW(igrid, tgrid, ygrid)]      =   amaxsave;
                                    searchWWE_J[inxW(igrid, tgrid, ygrid)]  =   searchval_J;
                                }
                                
                                if(iexcluded == 0){
                                    critereVF = max(critereVF,fabs(valueWWNE[inxW(igrid, tgrid, ygrid)] - (-valfnmax)));
                                    valueWWNE[inxW(igrid, tgrid, ygrid)]    =   -valfnmax;
                                    C_WWNE[inxW(igrid, tgrid, ygrid)]       =   params.C_val;
                                    saveWWNE[inxW(igrid, tgrid, ygrid)]     =   amaxsave;
                                    searchWWNE_J[inxW(igrid, tgrid, ygrid)] =   searchval_J;
                                }
                                
                            } // end iexcluded
                        } // for(ygrid)
                        
                    }//for(igrid=0;igrid<maxgrid;igrid++)
                }//for(tgrid=0;tgrid<maxtheta;tgrid++)
#if OMP == 1
            } // end pragma
#endif
            
            
            
            
            
#if indexPM == 0
            iter++;
            if(critereVF < crit_howard_off){howard_off = 1.0;}
            printf("Rule CNVG %d\t%20.15f\n",iter,critereVF);// getchar();
            
        }//while (critereVF > epsilonValue)
#endif
        
        
        
        
        
        
        /*********************************/
        /** PROJECTED METHOD SAVE FILES **/
        /*********************************/
        
        /******** SAVE IN A FILE FOR PROJECTED METHODS ******/
        // (see top of the EUP.h to see how) //
#if NOPOLICY == 1
        saveVFforPM(valueUiNE, C_UiNE, saveUiNE, valueUiE, C_UiE, saveUiE, valueUnNE, C_UnNE, saveUnNE, valueUnE, C_UnE, saveUnE, valueWWNE, C_WWNE, saveWWNE, valueWWE, C_WWE, saveWWE, valueJE, C_JE, saveJE, valueJNE, C_JNE, saveJNE);
#endif
        
#if SEA == 1
        saveVFforPM(valueUiNE, saveUiNE, valueUiE, saveUiE, valueUnNE, saveUnNE, valueUnE, saveUnE, valueWWNE, saveWWNE, valueWWE, saveWWE, valueJE, saveJE, valueJNE, saveJNE, valueJEi, saveJEi, valueJNEi, saveJNEi);
#endif
        
        
        
        // SAVE IN FILES BY OCCUPATION //
#if logprint == 1
        
        // initialize files
        FILE *saveoutfileJEval, *saveoutfileJEpol, *saveoutfileUnE, *saveoutfileUiE, *saveoutfileWWE, *saveoutfileJNEval, *saveoutfileJNEpol, *saveoutfileUnNE, *saveoutfileUiNE, *saveoutfileWWNE, *saveoutfileB, *saveoutfileR, *saveoutfileEhat;
        
        /**** ENTREPRENEUR (Ehat) ****/
        saveoutfileEhat=fopen(Ehatfile, "w");setbuf ( saveoutfileEhat, NULL );
        fprintf(saveoutfileEhat,"%5s\t%5s\t%20s\t%20s\t%20s\t%20s\t%20s\n","THETA", "ZPROD", "ASSET", "valueR", "C_R",  "saveR", "searchR");
        for(zgrid=0;zgrid<maxfirmtype;zgrid++){
            for(int tgrid = 0; tgrid < maxtheta ; tgrid++){
                for(igridCOH=0;igridCOH<maxgridCOH;igridCOH++){
                    fprintf(saveoutfileEhat,"%5d\t%5d\t%20.15f\t%20.15f\t%20.15f\t%20.15f\t%20.15f\n",tgrid, zgrid, gridCOH[igridCOH],valueEhat[inxRB(igridCOH,tgrid,zgrid)],C_Ehat[inxRB(igridCOH,tgrid,zgrid)],saveEhat[inxRB(igridCOH,tgrid,zgrid)],searchEhat[inxRB(igridCOH,tgrid,zgrid)]);
                }
            }
        }
        fclose(saveoutfileEhat);
        
        
        
        /**** ENTREPRENEURS ****/
        saveoutfileJNEpol=fopen(JNEpol, "w");saveoutfileJNEval=fopen(JNEval, "w");saveoutfileJEpol=fopen(JEpol, "w");saveoutfileJEval=fopen(JEval, "w");
        setbuf ( saveoutfileJNEpol, NULL ); setbuf ( saveoutfileJNEval, NULL ); setbuf ( saveoutfileJEpol, NULL ); setbuf ( saveoutfileJEval, NULL );
        
        fprintf(saveoutfileJEval,"%5s\t%5s\t%20s\t%20s\t%20s\t%20s\n","THETA", "Z", "ASSET", "valueJE", "C_JE", "collateralJE");
        fprintf(saveoutfileJNEval,"%5s\t%5s\t%5s\t%20s\t%20s\t%20s\t%20s\t%20s\t%20s\n","THETA", "Z", "I", "ASSET", "valueJNE", "C_JNE", "collateralJNE", "SloanJNE", "rateJNE");
        fprintf(saveoutfileJEpol,"%5s\t%5s\t%5s\t%20s\t%20s\t%20s\t%20s\n","THETA", "ZETA", "Z", "ASSET", "valueJEz", "saveJE", "searchJE");
        fprintf(saveoutfileJNEpol,"%5s\t%5s\t%5s\t%20s\t%20s\t%20s\t%20s\t%20s\n","THETA", "ZETA", "Z",  "ASSET","valueJNEz", "saveJNE","searchJNE","defaultJNE");
        for(int tgrid = 0; tgrid < maxtheta ; tgrid++){
            for(igrid=0;igrid<maxgrid;igrid++){
                for(int zgrid=0;zgrid<maxfirmtype;zgrid++) {
                    fprintf(saveoutfileJEval,"%5d\t%5d\t%20.15f\t%20.15f\t%20.15f\t%20.15f\n",tgrid, zgrid, grid[igrid],valueJE[inxE(igrid,tgrid,zgrid)],C_JE[inxE(igrid,tgrid,zgrid)],collateralJE[inxE(igrid,tgrid,zgrid)]);
                    fprintf(saveoutfileJNEval,"%5d\t%5d\t%5d\t%20.15f\t%20.15f\t%20.15f\t%20.15f\t%20.15f\t%20.15f\n",tgrid, zgrid,igrid, grid[igrid],valueJNE[inxE(igrid,tgrid, zgrid)],C_JNE[inxE(igrid,tgrid, zgrid)],collateralJNE[inxE(igrid,tgrid, zgrid)], sloanJNE[inxE(igrid,tgrid, zgrid)], rateJNE[inxE(igrid,tgrid, zgrid)]);
                    
                    for(int z=0;z<maxfirmtype;z++) {
                        fprintf(saveoutfileJEpol,"%5d\t%5d\t%5d\t%20.15f\t%20.15f\t%20.15f\t%20.15f\n",tgrid, zgrid, z, grid[igrid],valueJEz[inxEE(igrid,tgrid,zgrid,z)],saveJE[inxEE(igrid,tgrid,zgrid,z)], searchJE_W[inxEE(igrid,tgrid,zgrid,z)]);
                        fprintf(saveoutfileJNEpol,"%5d\t%5d\t%5d\t%20.15f\t%20.15f\t%20.15f\t%20.15f\t%20.15f\n",tgrid, zgrid, z, grid[igrid],valueJNEz[inxEE(igrid,tgrid,zgrid,z)],saveJNE[inxEE(igrid,tgrid,zgrid,z)], searchJNE_W[inxEE(igrid,tgrid,zgrid,z)], defaultJNE[inxEE(igrid,tgrid,zgrid,z)]);
                    }
                }
            }
        }
        fclose(saveoutfileJNEpol);fclose(saveoutfileJNEval);fclose(saveoutfileJEpol);fclose(saveoutfileJEval);
        
        
#if NOPOLICY == 1
        /**** ENTREPRENEUR (B) ****/
        saveoutfileB=fopen(Bfile, "w"); setbuf ( saveoutfileB, NULL );
        fprintf(saveoutfileB,"%5s\t%5s\t%20s\t%20s\t%20s\t%20s\t%20s\n","THETA", "ZPROD", "ASSET", "valueB", "C_B", "saveB", "searchB");
        for(zgrid=0;zgrid<maxfirmtype;zgrid++){
            for(int tgrid = 0; tgrid < maxtheta ; tgrid++){
                for(igridCOH=0;igridCOH<maxgridCOH;igridCOH++){
                    fprintf(saveoutfileB,"%5d\t%5d\t%20.15f\t%20.15f\t%20.15f\t%20.15f\t%20.15f\n",tgrid, zgrid, gridCOH[igridCOH],valueB[inxRB(igridCOH,tgrid,zgrid)],C_B[inxRB(igridCOH,tgrid,zgrid)],saveB[inxRB(igridCOH,tgrid,zgrid)],searchB[inxRB(igridCOH,tgrid,zgrid)]);
                }
            }
        }
        fclose(saveoutfileB);
        
        /**** ENTREPRENEUR (R) ****/
        saveoutfileR=fopen(Rfile, "w"); setbuf ( saveoutfileR, NULL );
        fprintf(saveoutfileR,"%5s\t%5s\t%20s\t%20s\t%20s\t%20s\t%20s\n","THETA", "ZPROD", "ASSET", "valueR", "C_R", "saveR", "searchR");
        for(zgrid=0;zgrid<maxfirmtype;zgrid++){
            for(int tgrid = 0; tgrid < maxtheta ; tgrid++){
                for(igridCOH=0;igridCOH<maxgridCOH;igridCOH++){
                    fprintf(saveoutfileR,"%5d\t%5d\t%20.15f\t%20.15f\t%20.15f\t%20.15f\t%20.15f\n", tgrid, zgrid, gridCOH[igridCOH],valueR[inxRB(igridCOH,tgrid,zgrid)],C_R[inxRB(igridCOH,tgrid,zgrid)],saveR[inxRB(igridCOH,tgrid,zgrid)],searchR[inxRB(igridCOH,tgrid,zgrid)]);
                }
            }
        }
        fclose(saveoutfileR);
#endif
        
        
#if SEA == 1
        /**** ENTREPRENEUR (B) ****/
        saveoutfileB=fopen(Bfile, "w"); setbuf ( saveoutfileB, NULL );
        fprintf(saveoutfileB,"%5s\t%5s\t%20s\t%20s\t%20s\t%20s\t%20s\t%20s\t%20s\t%20s\n","THETA", "ZPROD", "ASSET", "valueB", "C_B", "saveB", "searchB", "valueBi", "saveBi", "searchBi");
        for(zgrid=0;zgrid<maxfirmtype;zgrid++){
            for(int tgrid = 0; tgrid < maxtheta ; tgrid++){
                for(igridCOH=0;igridCOH<maxgridCOH;igridCOH++){
                    fprintf(saveoutfileB,"%5d\t%5d\t%20.15f\t%20.15f\t%20.15f\t%20.15f\t%20.15f\t%20.15f\t%20.15f\t%20.15f\n",tgrid, zgrid, gridCOH[igridCOH],valueB[inxRB(igridCOH,tgrid,zgrid)],C_B[inxRB(igridCOH,tgrid,zgrid)],saveB[inxRB(igridCOH,tgrid,zgrid)],searchB[inxRB(igridCOH,tgrid,zgrid)],valueBi[inxRB(igridCOH,tgrid,zgrid)],saveBi[inxRB(igridCOH,tgrid,zgrid)],searchBi[inxRB(igridCOH,tgrid,zgrid)]);
                }
            }
        }
        fclose(saveoutfileB);
        
        /**** ENTREPRENEUR (R) ****/
        saveoutfileR=fopen(Rfile, "w"); setbuf ( saveoutfileR, NULL );
        fprintf(saveoutfileR,"%5s\t%5s\t%20s\t%20s\t%20s\t%20s\t%20s\t%20s\t%20s\t%20s\n","THETA", "ZPROD", "ASSET", "valueR", "C_R",  "saveR", "searchR", "valueRi", "saveRi", "searchRi");
        for(zgrid=0;zgrid<maxfirmtype;zgrid++){
            for(int tgrid = 0; tgrid < maxtheta ; tgrid++){
                for(igridCOH=0;igridCOH<maxgridCOH;igridCOH++){
                    fprintf(saveoutfileR,"%5d\t%5d\t%20.15f\t%20.15f\t%20.15f\t%20.15f\t%20.15f\t%20.15f\t%20.15f\t%20.15f\n", tgrid, zgrid, gridCOH[igridCOH],valueR[inxRB(igridCOH,tgrid,zgrid)],C_R[inxRB(igridCOH,tgrid,zgrid)],saveR[inxRB(igridCOH,tgrid,zgrid)],searchR[inxRB(igridCOH,tgrid,zgrid)],valueRi[inxRB(igridCOH,tgrid,zgrid)],saveRi[inxRB(igridCOH,tgrid,zgrid)],searchRi[inxRB(igridCOH,tgrid,zgrid)]);
                }
            }
        }
        fclose(saveoutfileR);
        
        
        saveoutfileJNEpol=fopen(JNEipol, "w");saveoutfileJNEval=fopen(JNEival, "w");saveoutfileJEpol=fopen(JEipol, "w");saveoutfileJEval=fopen(JEival, "w");
        setbuf ( saveoutfileJNEpol, NULL );setbuf ( saveoutfileJNEval, NULL );setbuf ( saveoutfileJEpol, NULL );setbuf ( saveoutfileJEval, NULL );
        
        fprintf(saveoutfileJEval,"%5s\t%5s\t%20s\t%20s\t%20s\t%20s\n","THETA", "Z", "ASSET", "valueJE", "C_JE", "collateralJE");
        fprintf(saveoutfileJNEval,"%5s\t%5s\t%20s\t%20s\t%20s\t%20s\t%20s\t%20s\n","THETA", "Z", "ASSET", "valueJNE", "C_JNE", "collateralJNE", "SloanJNE", "rateJNE");
        fprintf(saveoutfileJEpol,"%5s\t%5s\t%5s\t%20s\t%20s\t%20s\t%20s\n","THETA", "ZETA", "Z", "ASSET", "valueJEz", "saveJE", "searchJE");
        fprintf(saveoutfileJNEpol,"%5s\t%5s\t%5s\t%20s\t%20s\t%20s\t%20s\t%20s\n","THETA", "ZETA", "Z",  "ASSET","valueJNEz", "saveJNE","searchJNE","defaultJNE");
        for(int zgrid=0;zgrid<maxfirmtype;zgrid++) {
            for(int tgrid = 0; tgrid < maxtheta ; tgrid++)
            {
                for(igrid=0;igrid<maxgrid;igrid++)
                {
                    for(int zgrid=0;zgrid<maxfirmtype;zgrid++) {
                        fprintf(saveoutfileJEval,"%5d\t%5d\t%20.15f\t%20.15f\t%20.15f\t%20.15f\n",tgrid, zgrid, grid[igrid],valueJEi[inxE(igrid,tgrid,zgrid)], C_JEi[inxE(igrid,tgrid,zgrid)], collateralJEi[inxE(igrid,tgrid,zgrid)]);
                        fprintf(saveoutfileJNEval,"%5d\t%5d\t%20.15f\t%20.15f\t%20.15f\t%20.15f\t%20.15f\t%20.15f\n",tgrid, zgrid, grid[igrid],valueJNEi[inxE(igrid,tgrid, zgrid)], C_JNEi[inxE(igrid,tgrid, zgrid)], collateralJNEi[inxE(igrid,tgrid, zgrid)], sloanJNEi[inxE(igrid,tgrid, zgrid)], rateJNEi[inxE(igrid,tgrid, zgrid)]);
                        
                        for(int z=0;z<maxfirmtype;z++) {
                            fprintf(saveoutfileJEpol,"%5d\t%5d\t%5d\t%20.15f\t%20.15f\t%20.15f\t%20.15f\n",tgrid, zgrid, z, grid[igrid],valueJEiz[inxEE(igrid,tgrid,zgrid,z)],saveJEi[inxEE(igrid,tgrid,zgrid,z)], searchJEi_W[inxEE(igrid,tgrid,zgrid,z)]);
                            fprintf(saveoutfileJNEpol,"%5d\t%5d\t%5d\t%20.15f\t%20.15f\t%20.15f\t%20.15f\t%20.15f\n",tgrid, zgrid, z, grid[igrid],valueJNEiz[inxEE(igrid,tgrid,zgrid,z)],saveJNEi[inxEE(igrid,tgrid,zgrid,z)], searchJNEi_W[inxEE(igrid,tgrid,zgrid,z)], defaultJNEi[inxEE(igrid,tgrid,zgrid,z)]);
                        }
                    }
                }
            }
        }
        
        
        fclose(saveoutfileJNEpol);fclose(saveoutfileJNEval);fclose(saveoutfileJEpol);fclose(saveoutfileJEval);
#endif
        
        
        /**** UNEMPLOYED ****/
        saveoutfileUnE=fopen(UnE, "w");saveoutfileUiE=fopen(UiE, "w");saveoutfileUnNE=fopen(UnNE, "w");saveoutfileUiNE=fopen(UiNE, "w");
        setbuf ( saveoutfileUnE, NULL );setbuf ( saveoutfileUiE, NULL );setbuf ( saveoutfileUnNE, NULL );setbuf ( saveoutfileUiNE, NULL );
        
        fprintf(saveoutfileUnE,"%5s\t%20s\t%20s\t%20s\t%20s\t%20s\t%20s\n","THETA", "ASSET", "valueUnE", "C_UnE", "saveUnE", "searchUnE_J", "searchUnE_W");
        fprintf(saveoutfileUiE,"%5s\t%20s\t%20s\t%20s\t%20s\t%20s\t%20s\n","THETA", "ASSET", "valueUiE", "C_UiE", "saveUiE", "searchUiE_J", "searchUiE_W");
        fprintf(saveoutfileUnNE,"%5s\t%20s\t%20s\t%20s\t%20s\t%20s\t%20s\n","THETA", "ASSET", "valueUnNE", "C_UnNE" "saveUnNE", "searchUnNE_J", "searchUnNE_W");
        fprintf(saveoutfileUiNE,"%5s\t%20s\t%20s\t%20s\t%20s\t%20s\t%20s\n","THETA", "ASSET", "valueUiNE", "C_UiNE", "saveUiNE", "searchUiNE_J", "searchUiNE_W");
        
        for(int tgrid = 0; tgrid < maxtheta ; tgrid++){
            for(igrid=0;igrid<maxgrid;igrid++){
                fprintf(saveoutfileUnE,"%5d\t%20.15f\t%20.15f\t%20.15f\t%20.15f\t%20.15f\t%20.15f\n",tgrid, grid[igrid],valueUnE[inx(igrid,tgrid)],C_UnE[inx(igrid,tgrid)],saveUnE[inx(igrid,tgrid)],searchUnE_J[inx(igrid,tgrid)],searchUnE_W[inx(igrid,tgrid)]);
                fprintf(saveoutfileUiE,"%5d\t%20.15f\t%20.15f\t%20.15f\t%20.15f\t%20.15f\t%20.15f\n",tgrid, grid[igrid],valueUiE[inx(igrid,tgrid)],C_UiE[inx(igrid,tgrid)],saveUiE[inx(igrid,tgrid)],searchUiE_J[inx(igrid,tgrid)],searchUiE_W[inx(igrid,tgrid)]);
                fprintf(saveoutfileUnNE,"%5d\t%20.15f\t%20.15f\t%20.15f\t%20.15f\t%20.15f\t%20.15f\n",tgrid, grid[igrid],valueUnNE[inx(igrid,tgrid)],C_UnNE[inx(igrid,tgrid)],saveUnNE[inx(igrid,tgrid)],searchUnNE_J[inx(igrid,tgrid)],searchUnNE_W[inx(igrid,tgrid)]);
                fprintf(saveoutfileUiNE,"%5d\t%20.15f\t%20.15f\t%20.15f\t%20.15f\t%20.15f\t%20.15f\n",tgrid, grid[igrid],valueUiNE[inx(igrid,tgrid)],C_UiNE[inx(igrid,tgrid)],saveUiNE[inx(igrid,tgrid)],searchUiNE_J[inx(igrid,tgrid)],searchUiNE_W[inx(igrid,tgrid)]);
            }
        }
        fclose(saveoutfileUiE);fclose(saveoutfileUiNE);fclose(saveoutfileUnE);fclose(saveoutfileUnNE);
        
        
        /**** WORKER ****/
        saveoutfileWWE=fopen(WWE, "w");saveoutfileWWNE=fopen(WWNE, "w");
        setbuf ( saveoutfileWWE, NULL );setbuf ( saveoutfileWWNE, NULL );
        fprintf(saveoutfileWWNE,"%5s\t%5s\t%20s\t%20s\t%20s\t%20s\t%20s\n","THETA", "PROD", "ASSET", "valueWWNE", "C_WWNE", "saveWWNE", "searchWWNE_J");
        fprintf(saveoutfileWWE,"%5s\t%5s\t%20s\t%20s\t%20s\t%20s\t%20s\n","THETA", "PROD", "ASSET", "valueWWE", "C_WWE", "saveWWE", "searchWWE_J");
        for(int tgrid = 0; tgrid < maxtheta ; tgrid++){
            for(int ygrid = 0; ygrid < maxyprod; ygrid++){
                for(igrid=0;igrid<maxgrid;igrid++){
                    fprintf(saveoutfileWWE,"%5d\t%5d\t%20.15f\t%20.15f\t%20.15f\t%20.15f\t%20.15f\n",tgrid, ygrid, grid[igrid],valueWWE[inxW(igrid,tgrid, ygrid)],C_WWE[inxW(igrid,tgrid, ygrid)],saveWWE[inxW(igrid,tgrid, ygrid)],searchWWE_J[inxW(igrid,tgrid, ygrid)]);
                    fprintf(saveoutfileWWNE,"%5d\t%5d\t%20.15f\t%20.15f\t%20.15f\t%20.15f\t%20.15f\n",tgrid, ygrid, grid[igrid],valueWWNE[inxW(igrid,tgrid, ygrid)],C_WWNE[inxW(igrid,tgrid, ygrid)],saveWWNE[inxW(igrid,tgrid, ygrid)],searchWWNE_J[inxW(igrid,tgrid, ygrid)]);
                }
            }
        }
        fclose(saveoutfileWWE);fclose(saveoutfileWWNE);
        
        
        /**** UNEMPLOYED SEARCH EFFORTS ****/
#if writeSESW
        FILE *SWbenchUiE, *SWbenchUiNE, *SEbenchUiE, *SEbenchUiNE;
#if writebenchSESW == 1
        SEbenchUiE=fopen(searchSEUiEbench, "w");SWbenchUiE=fopen(searchSWUiEbench, "w");SEbenchUiNE=fopen(searchSEUiNEbench, "w");SWbenchUiNE=fopen(searchSWUiNEbench, "w");
#endif
#if writereformSESW == 1
        SEbenchUiE=fopen(searchSEUiEreform, "w");SWbenchUiE=fopen(searchSWUiEreform, "w");SEbenchUiNE=fopen(searchSEUiNEreform, "w");SWbenchUiNE=fopen(searchSWUiNEreform, "w");
#endif
        setbuf ( SEbenchUiE, NULL );setbuf ( SWbenchUiE, NULL );setbuf ( SEbenchUiNE, NULL );setbuf ( SWbenchUiNE, NULL );
        for(igrid=0;igrid<maxgrid;igrid++){
            for(int tgrid = 0; tgrid < maxtheta ; tgrid++){
                fprintf(SEbenchUiE,"%20.15f\n",searchUiE_J[inx(igrid,tgrid)]);
                fprintf(SWbenchUiE,"%20.15f\n",searchUiE_W[inx(igrid,tgrid)]);
                fprintf(SEbenchUiNE,"%20.15f\n",searchUiNE_J[inx(igrid,tgrid)]);
                fprintf(SWbenchUiNE,"%20.15f\n",searchUiNE_W[inx(igrid,tgrid)]);
            }
        }
        fclose(SEbenchUiE);fclose(SWbenchUiE);fclose(SEbenchUiNE);fclose(SWbenchUiNE);
#endif
#endif // logprint
        
        
        
        /*******************/
        /** FREE POINTERS **/
        /*******************/
#if OPT_compute_C_flow == 1
        free(C_expE_UnE); free(C_expE_UiE); free(C_expE_WWE); free(C_expE_WWNE_JNE); free(C_expE_WWE_JE); free(C_expE_UnNE_JNE); free(C_expE_UnE_JE);
        free(C_expU_UnNE_JNE); free(C_expU_UiNE_JNE); free(C_expU_UnE_JE); free(C_expU_UiE_JE); free(C_expU_UnNE); free(C_expU_UnE); free(C_expU_UiNE); free(C_expU_UiE); free(C_expU_WWNE); free(C_expU_WWE); free(C_expU_WWNE_JNE); free(C_expU_WWE_JE);
        free(C_expW_WWE_JE); free(C_expW_WWNE_JNE); free(C_expW_WWE); free(C_expW_WWNE);
        free(C_B); free(C_R); free(C_Ehat);
#endif
        free(valueJNEtemp); free(saveJNEtemp); free(searchJNEtemp); free(defaultJNEtemp); free(valueJEtemp); free(saveJEtemp); free(searchJEtemp);
        free(valueB); free(valueR);free(valueEhat);
        free(valueJEz); free(valueJNEz);
        free(searchB);free(searchR);free(searchEhat);free(saveB);free(saveR);free(saveEhat);
        
        free(expE_UnE); free(expE_UiE); free(expE_WWE); free(expE_WWNE_JNE); free(expE_WWE_JE); free(expE_UnNE_JNE); free(expE_UnE_JE);
        free(expU_UnNE_JNE); free(expU_UiNE_JNE); free(expU_UnE_JE); free(expU_UiE_JE); free(expU_UnNE); free(expU_UnE); free(expU_UiNE);
        free(expU_UiE); free(expU_WWNE); free(expU_WWE); free(expU_WWNE_JNE); free(expU_WWE_JE);
        free(expW_WWE_JE); free(expW_WWNE_JNE); free(expW_WWE); free(expW_WWNE);
        
#if SEA == 1
#if OPT_compute_C_flow == 1
        free(C_expE_WWNE_JNEi); free(C_expE_WWE_JEi); free(C_expE_UiE_JEi); free(C_expE_UiNE_JNEi);
        free(C_expU_UiNE_JNEi); free(C_expU_UiE_JEi); free(C_expU_WWNE_JNEi); free(C_expU_WWE_JEi);
#endif
        free(valueBi); free(valueRi);free(valueEhati);
        free(C_Bi); free(C_Ri); free(C_Ehati);
        free(searchBi);free(searchRi);free(searchEhati);free(saveBi);free(saveRi);free(saveEhati);
        free(valueJEiz); free(valueJNEiz);
        free(expE_WWNE_JNEi); free(expE_WWE_JEi); free(expE_UiE_JEi); free(expE_UiNE_JNEi);
        free(expU_UiNE_JNEi); free(expU_UiE_JEi); free(expU_WWNE_JNEi); free(expU_WWE_JEi);
#endif
        
    }
